{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4dc1bf63",
   "metadata": {},
   "source": [
    "# 11 Non-Compliance and Instruments\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a6d9bc21",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:12.624058Z",
     "start_time": "2024-01-23T12:22:09.270731Z"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "from toolz.curried import *\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from scipy.special import expit\n",
    "\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "import seaborn as sns\n",
    "from matplotlib import pyplot as plt\n",
    "import matplotlib\n",
    "\n",
    "\n",
    "from cycler import cycler\n",
    "\n",
    "color=['0.0', '0.4', '0.8']\n",
    "default_cycler = (cycler(color=color))\n",
    "linestyle=['-', '--', ':', '-.']\n",
    "marker=['o', 'v', 'd', 'p']\n",
    "\n",
    "plt.rc('axes', prop_cycle=default_cycler)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f0a112e",
   "metadata": {},
   "source": [
    "## Non-Compliance\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4a0845f8",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:13.194109Z",
     "start_time": "2024-01-23T12:22:12.625731Z"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.50.0 (0)\n",
       " -->\n",
       "<!-- Pages: 1 -->\n",
       "<svg width=\"242pt\" height=\"98pt\"\n",
       " viewBox=\"0.00 0.00 242.00 98.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 94)\">\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-94 238,-94 238,4 -4,4\"/>\n",
       "<!-- U -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>U</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">U</text>\n",
       "</g>\n",
       "<!-- T -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>T</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"117\" cy=\"-64\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"117\" y=\"-60.3\" font-family=\"Times,serif\" font-size=\"14.00\">T</text>\n",
       "</g>\n",
       "<!-- U&#45;&gt;T -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>U&#45;&gt;T</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M48.85,-28.9C60.03,-34.74 74,-42.04 86.29,-48.47\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"84.81,-51.65 95.29,-53.18 88.05,-45.44 84.81,-51.65\"/>\n",
       "</g>\n",
       "<!-- Y -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>Y</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"207\" cy=\"-41\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"207\" y=\"-37.3\" font-family=\"Times,serif\" font-size=\"14.00\">Y</text>\n",
       "</g>\n",
       "<!-- U&#45;&gt;Y -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>U&#45;&gt;Y</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M53.88,-21.34C84.55,-25.31 135.83,-31.93 170.27,-36.38\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"170.01,-39.88 180.37,-37.69 170.91,-32.94 170.01,-39.88\"/>\n",
       "</g>\n",
       "<!-- T&#45;&gt;Y -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>T&#45;&gt;Y</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M142.51,-57.59C151.57,-55.23 162.01,-52.5 171.78,-49.95\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"172.83,-53.29 181.62,-47.37 171.06,-46.51 172.83,-53.29\"/>\n",
       "</g>\n",
       "<!-- Z -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>Z</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-72\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-68.3\" font-family=\"Times,serif\" font-size=\"14.00\">Z</text>\n",
       "</g>\n",
       "<!-- Z&#45;&gt;T -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>Z&#45;&gt;T</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M53.93,-69.64C62.09,-68.9 71.26,-68.07 79.99,-67.27\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"80.31,-70.76 89.95,-66.37 79.68,-63.79 80.31,-70.76\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.graphs.Digraph at 0x7f84bb4da950>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from graphviz import Digraph\n",
    "\n",
    "gr = Digraph(format=\"png\", graph_attr={\"rankdir\":\"LR\"})\n",
    "\n",
    "gr.edge(\"U\", \"T\")\n",
    "gr.edge(\"U\", \"Y\")\n",
    "gr.edge(\"Z\", \"T\")\n",
    "gr.edge(\"T\", \"Y\")\n",
    "\n",
    "gr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "46e508be",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:13.225362Z",
     "start_time": "2024-01-23T12:22:13.195735Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>age</th>\n",
       "      <th>income</th>\n",
       "      <th>credit_score</th>\n",
       "      <th>prime_elegible</th>\n",
       "      <th>prime_card</th>\n",
       "      <th>pv</th>\n",
       "      <th>tau</th>\n",
       "      <th>categ</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>37.7</td>\n",
       "      <td>9687.0</td>\n",
       "      <td>822.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>4913.79</td>\n",
       "      <td>700.0</td>\n",
       "      <td>complier</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>46.0</td>\n",
       "      <td>13731.0</td>\n",
       "      <td>190.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>5637.66</td>\n",
       "      <td>200.0</td>\n",
       "      <td>never-taker</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>43.1</td>\n",
       "      <td>2839.0</td>\n",
       "      <td>214.0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2410.45</td>\n",
       "      <td>700.0</td>\n",
       "      <td>complier</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>36.0</td>\n",
       "      <td>1206.0</td>\n",
       "      <td>318.0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1363.06</td>\n",
       "      <td>700.0</td>\n",
       "      <td>complier</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>39.7</td>\n",
       "      <td>4095.0</td>\n",
       "      <td>430.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2189.80</td>\n",
       "      <td>700.0</td>\n",
       "      <td>complier</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    age   income  credit_score  prime_elegible  prime_card       pv    tau  \\\n",
       "0  37.7   9687.0         822.0               0           0  4913.79  700.0   \n",
       "1  46.0  13731.0         190.0               0           0  5637.66  200.0   \n",
       "2  43.1   2839.0         214.0               1           1  2410.45  700.0   \n",
       "3  36.0   1206.0         318.0               1           1  1363.06  700.0   \n",
       "4  39.7   4095.0         430.0               0           0  2189.80  700.0   \n",
       "\n",
       "         categ  \n",
       "0     complier  \n",
       "1  never-taker  \n",
       "2     complier  \n",
       "3     complier  \n",
       "4     complier  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "df = pd.read_csv(\"./data/prime_card.csv\")\n",
    "\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f266f971",
   "metadata": {},
   "source": [
    "## Extending Potential Outcomes\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6b4146eb",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:13.249112Z",
     "start_time": "2024-01-23T12:22:13.227557Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "         <td></td>           <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>      <td> 2498.3618</td> <td>   24.327</td> <td>  102.701</td> <td> 0.000</td> <td> 2450.677</td> <td> 2546.047</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>prime_elegible</th> <td>  321.3880</td> <td>   34.321</td> <td>    9.364</td> <td> 0.000</td> <td>  254.113</td> <td>  388.663</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = smf.ols(\"pv~prime_elegible\", data=df).fit()\n",
    "m.summary().tables[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "bf978a12",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:13.254232Z",
     "start_time": "2024-01-23T12:22:13.250451Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "413.45"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df[\"tau\"].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "179a99c7",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:13.268017Z",
     "start_time": "2024-01-23T12:22:13.255653Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "       <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>  <td> 2534.4947</td> <td>   19.239</td> <td>  131.740</td> <td> 0.000</td> <td> 2496.783</td> <td> 2572.206</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>prime_card</th> <td>  588.1388</td> <td>   41.676</td> <td>   14.112</td> <td> 0.000</td> <td>  506.446</td> <td>  669.831</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = smf.ols(\"pv~prime_card\", data=df).fit()\n",
    "m.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb010ef5",
   "metadata": {},
   "source": [
    "## Instrument Identification Assumptions\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c90cca32",
   "metadata": {},
   "source": [
    "## First Stage\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "471dfd47",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:13.281666Z",
     "start_time": "2024-01-23T12:22:13.269945Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "         <td></td>           <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>      <td> 6.729e-15</td> <td>    0.005</td> <td> 1.35e-12</td> <td> 1.000</td> <td>   -0.010</td> <td>    0.010</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>prime_elegible</th> <td>    0.4242</td> <td>    0.007</td> <td>   60.536</td> <td> 0.000</td> <td>    0.410</td> <td>    0.438</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "first_stage = smf.ols(\"prime_card ~ prime_elegible\", data=df).fit()\n",
    "first_stage.summary().tables[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ba1daf77",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:13.292213Z",
     "start_time": "2024-01-23T12:22:13.283317Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "categ\n",
       "complier       0.4269\n",
       "never-taker    0.5731\n",
       "dtype: float64"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.groupby(\"categ\").size()/len(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3672e47d",
   "metadata": {},
   "source": [
    "## Reduced Form\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d4e10cdf",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:13.306388Z",
     "start_time": "2024-01-23T12:22:13.293975Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "         <td></td>           <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>      <td> 2498.3618</td> <td>   24.327</td> <td>  102.701</td> <td> 0.000</td> <td> 2450.677</td> <td> 2546.047</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>prime_elegible</th> <td>  321.3880</td> <td>   34.321</td> <td>    9.364</td> <td> 0.000</td> <td>  254.113</td> <td>  388.663</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "red_form = smf.ols(\"pv ~ prime_elegible\", data=df).fit()\n",
    "red_form.summary().tables[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f583cf79",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:13.313996Z",
     "start_time": "2024-01-23T12:22:13.310505Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "757.6973795343938"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "late = (red_form.params[\"prime_elegible\"] /\n",
    "        first_stage.params[\"prime_elegible\"])\n",
    "late"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "c1d57b18",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:13.322436Z",
     "start_time": "2024-01-23T12:22:13.315622Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "categ\n",
       "complier       700.0\n",
       "never-taker    200.0\n",
       "Name: tau, dtype: float64"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.groupby(\"categ\")[\"tau\"].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d63f18a1",
   "metadata": {},
   "source": [
    "## Two Stage Least Squares\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f44cd90b",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:13.478692Z",
     "start_time": "2024-01-23T12:22:13.323966Z"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.50.0 (0)\n",
       " -->\n",
       "<!-- Pages: 1 -->\n",
       "<svg width=\"152pt\" height=\"98pt\"\n",
       " viewBox=\"0.00 0.00 152.00 98.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 94)\">\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-94 148,-94 148,4 -4,4\"/>\n",
       "<!-- U -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>U</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">U</text>\n",
       "</g>\n",
       "<!-- T -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>T</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"117\" cy=\"-72\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"117\" y=\"-68.3\" font-family=\"Times,serif\" font-size=\"14.00\">T</text>\n",
       "</g>\n",
       "<!-- U&#45;&gt;T -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>U&#45;&gt;T</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M47.53,-29.98C59.38,-37.26 74.7,-46.66 87.84,-54.72\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"86.17,-57.8 96.52,-60.05 89.83,-51.84 86.17,-57.8\"/>\n",
       "</g>\n",
       "<!-- Y -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>Y</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"117\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"117\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">Y</text>\n",
       "</g>\n",
       "<!-- U&#45;&gt;Y -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>U&#45;&gt;Y</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M54.4,-18C62.39,-18 71.31,-18 79.82,-18\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"79.92,-21.5 89.92,-18 79.92,-14.5 79.92,-21.5\"/>\n",
       "</g>\n",
       "<!-- Z -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>Z</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-72\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-68.3\" font-family=\"Times,serif\" font-size=\"14.00\">Z</text>\n",
       "</g>\n",
       "<!-- Z&#45;&gt;T -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>Z&#45;&gt;T</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M54.4,-72C62.39,-72 71.31,-72 79.82,-72\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"79.92,-75.5 89.92,-72 79.92,-68.5 79.92,-75.5\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.graphs.Digraph at 0x7f84882c5c90>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gr = Digraph(format=\"png\", graph_attr={\"rankdir\":\"LR\"})\n",
    "\n",
    "gr.edge(\"U\", \"T\")\n",
    "gr.edge(\"U\", \"Y\")\n",
    "gr.edge(\"Z\", \"T\")\n",
    "gr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "9607bf02",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:13.498654Z",
     "start_time": "2024-01-23T12:22:13.480857Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "       <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>  <td> 2498.3618</td> <td>   24.327</td> <td>  102.701</td> <td> 0.000</td> <td> 2450.677</td> <td> 2546.047</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>prime_card</th> <td>  757.6974</td> <td>   80.914</td> <td>    9.364</td> <td> 0.000</td> <td>  599.091</td> <td>  916.304</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "iv_regr = smf.ols(\n",
    "    \"pv ~ prime_card\",\n",
    "    data=df.assign(prime_card=first_stage.fittedvalues)).fit()\n",
    "\n",
    "iv_regr.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c851c3c",
   "metadata": {},
   "source": [
    "## Standard Error\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "e99ccfde",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:13.509300Z",
     "start_time": "2024-01-23T12:22:13.500116Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SE IV: 80.52861026141942\n",
      "95% CI: [596.6401590115549, 918.7546000572327]\n"
     ]
    }
   ],
   "source": [
    "Z = df[\"prime_elegible\"]\n",
    "T = df[\"prime_card\"]\n",
    "n = len(df)\n",
    "\n",
    "# not the same as iv_regr.resid!\n",
    "e_iv = df[\"pv\"] - iv_regr.predict(df)\n",
    "compliance = np.cov(T, Z)[0, 1]/Z.var()\n",
    "\n",
    "se = np.std(e_iv)/(compliance*np.std(Z)*np.sqrt(n))\n",
    "\n",
    "print(\"SE IV:\", se)\n",
    "print(\"95% CI:\", [late - 2*se, late + 2*se])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "55102be8",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:13.659991Z",
     "start_time": "2024-01-23T12:22:13.511367Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Parameter Estimates</caption>\n",
       "<tr>\n",
       "       <td></td>      <th>Parameter</th> <th>Std. Err.</th> <th>T-stat</th> <th>P-value</th> <th>Lower CI</th> <th>Upper CI</th>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>   <td>2498.4</td>    <td>24.211</td>   <td>103.19</td> <td>0.0000</td>   <td>2450.9</td>   <td>2545.8</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>prime_card</th>  <td>757.70</td>    <td>80.529</td>   <td>9.4090</td> <td>0.0000</td>   <td>599.86</td>   <td>915.53</td> \n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from linearmodels import IV2SLS\n",
    "\n",
    "formula = 'pv ~ 1 + [prime_card ~ prime_elegible]'\n",
    "iv_model = IV2SLS.from_formula(formula, df).fit(cov_type=\"unadjusted\")\n",
    "\n",
    "iv_model.summary.tables[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "e7f8cb91",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:14.017449Z",
     "start_time": "2024-01-23T12:22:13.662388Z"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f84bb78d8d0>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAz4AAAFNCAYAAADB4rqeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABpAElEQVR4nO3dd3gU5frG8e+TEAi9g0Bo0kuQEhWwrTQbAgooqAiKothQ5CjYPTaOB7AcFQUPRiyAtJ8UFRApAiIdQbqgEHpAeg/v74/d7EmogZTZJPfnuvbKvjOzu/cMYSbPzjvvmHMOERERERGRrCzM6wAiIiIiIiLpTYWPiIiIiIhkeSp8REREREQky1PhIyIiIiIiWZ4KHxERERERyfJU+IiIiIiISJanwkckCzGzWDN73escIiLZhZk9Z2afpuH7TTezB9Lq/dKDmVUwM2dmOdLhvQ+Y2aVp/b5pycx+NzOf1znkwqnwEc+YWQ0z+8nM9prZOjO7Lcm8xJ3qgSSPF5PMv8vMtprZhqQ7HzOrZGZzzCz8PJ9dysz+G3iP/Wa2ysxeNbO8gfnOzCqn/Vp7x8y6mNmsFCzjzOyOQPuaJNv/4Bn+TcoFDtJHTpk+PmPWSkTEz8z+NLPDgX3QtsAXQfnS+3Odc2865zKkUDGzQmY2JLB++81sjZk9mxGfnVbOtw7OuXzOufVp+Hmn/j3xp5n1voDXn/aFonOulnNuelpllIyjwkc8EfiW6FtgAlAE6AZ8aWZVT1m0UGAnmM8591qS1/YF6gOPAx8kWf59oKdzLuEcn10E+AXIDTRyzuUHmgOFgEppsHpJPyvNvw1LZ52B3YGfOOd+Ttz+QK3AMkn/TTYGpj2WZFo+59ytHmQXEbk1sL+qC9QD+ngZJh2OAe8A+YAaQEGgFfBHGn9GevNqHQoFfjfaAS+aWfMM+EwJMSp8xCvVgdLAO865BOfcT8BsoFMKXlsU2Oyc2wr8CFwKYGbtAtPnnuf1PYH9wD3OuT8BnHObnHM9nHO/ne/DzewWM1tsZvvMbJOZvZJkXuI3S13NbCPwk5mFmdkLZvaXme0ws6FmVjCwvM/M4k55/z/NrFng+Stm9k3gNfsDp9djkixbz8wWBeaNACLPl/8c61UeuA5/EXqDmZW82PcSEfGSc24bMAl/AQSAmTUM9AjYY2ZLT+ktUNHMZgT2pVPM7AMz+zIwLyX76cRlTzsGBKbfb2YrzexvM5sU2N8mvlfzQK+DvWb2AWDnWLXLga+dc387504651Y550Ylea/3AselfWa20MyuSTLvFTMbaWZfBtZzmZlVNbM+gWPTJjNrkWT56Wb2lpnNC2T7NvDF4WnMrKD9rxfFZjN73c7e8+J86+DMrLKZlbbkPQkOmZlLstxZt+m5OOcWAL+T/HdjpPnPQO01s5lmViswvRtwN/CMJenNcMq/fy4ze9fMtgQe75pZrpRkkYynwke8cqYduwG1T5n2l5nFmdlnZlYsMG0nUNTMovCfqfnd/N0ZXiBl3+41A8Y4505eZPaDwL34zxDdAnQ3szanLHMd/m+zbgC6BB7X4y/S8pH8LNX5tAKGBz5vXOJrzSwn8H/AF/jPmo0E2l7guiR1L7DAOTcaWIl/Zy8ikukEjg83AesC7TLAROB1/PvLXsBoMyseeMnXwEKgGPAagbPeqRA8BgSOD88BtwPFgZ+BYYFcxYDR+I9fxfCf+bjqHO87F3jDzO4zsypnmD8f/x/0RQLrNNLMkn4hdiv+Y0ZhYDH+4jAMKAP8E/jklPe7F7gf/xeVJ/D3qjiTzwPzK+M/09YCOFv3v/OtAwDOuS1JexIAY/EfCznXNj0fM2uI/2+NdUkmfw9UAUoAi4CvAhkGBZ6/fY7eDM8DDfFv98uAK/D/e0oocs7poUeGP4AIYD3wTOB5C+AYMCkwPx8QA+QASgKjEucF5jfFv/OcgX9nMwDoCviAafh35rXP8tlrgYfPk88BlVO4Lu/iP3MFUCHw2kuTzJ8KPJKkXQ04Hlg3HxB3yvv9CTQLPH8F+DHJvJrA4cDza4EtgCWZPwd4/Sw5uwCzzrEea4EnA8/7AEtPmZ+4bjlOmT4dOATsSfJ4zevfMT300CN7PQL7zgP4z+i7wL63UGDes8AXpyw/CX+BUw7/H+15k8z7Gvgy8Dwl++nEZc90DPge6JqkHRbYZ5bHX1jMTTLPgDjggbOsY278f/AvDBxH1gE3nWOb/A1cliTnlCTzbg1sr/BAO38ge+I2mw70TbJ8TfzH6fCkxwP8x+ijQO4ky3YEpl3MOnCG42/g329h4meca5ue4fMSs+4BDgee9yPJsfOU5QsFlikYaMdyynH1lH//P4Cbk8y7AfjT6/8Pepz5oTM+4gnn3HGgDf4zJtuAp4Fv8O/wcc4dcM4tcM6dcM5tBx4DWphZgcD8qc65hs6564CT+IukWPzfZHXB/43d2UbZ2QWUutjsZnalmU0zs51mthd4GP83dUltSvK8NPBXkvZf/O9gkRLbkjw/BESav994afxd+1yS+X9xEczsKqAigW/T8B/0o82sbgrf4gnnXKEkjxfP/xIRkTTXxvmv2/Th71KduG8uD7QPdHPbY2Z7gKvxHwtKA3875w4meZ+L2pcmkfQYUB54L8nn7sZf4JQJfHZw2cD+POlrk3HOHXb+wRQa4O/2/Q3+szpFAMzs6UD3r72BzypI8uPT9iTPDwPx7n/XxB4O/Ew6IETSLH/h/6Ly1ONd+cD0rUnW8RP8Z08ueB1OZWY3AT3w/9smZjzXNj2bYoF164X/9yMi8P7hZtbXzP4ws334i5rE5VPiTMf40il8rWQwFT7iGefcb86565xzRZ1zN+DvBjbvbIsHfibrImdmhr/r1xP4d1Lhzrm/8J/ur3OW9/oRuM3MLvb3/2v8Xc7KOucKAh+fmitJXvCflUna9zjx28Xt+LvN5UmyPuH4T9unxFagTGAbJH3vi9EZ/zosMbNtwK+B6fde5PuJiHjGOTcD/5dh/QKTNuE/45P0C5q8zrm++PelhS0wqmdA0n3pxeynkx4DNgEPnfLZuZ1zcwKfXTbJe1vS9nnWcR/wJpAXqBi4nudZ4A6gsHOuELCXc18zdD5Js5TDf4Ym/pRlNuE/41MsyfoVcM7V4jxOXYdT55tZNfzd6O5wziUtws61Tc/1eQnOuf7AEeCRwOS7gNb4u8EXxH+GCP633ZL+W57JmY7xW87zGvGICh/xjJnVMbNIM8tjZr3wf/MWG5h3pZlVM//AAEXx9yue7pzbe8rbPAAsds4twX8mJ7eZ1cR/Pc3ZhsMcABQAPk+8GNLMypjZADM7W7GUVH5gt3PuiJldgX+neS7DgKfMf/FsPvw7+RHOuRPAGvxncG4xswj8/YJTelHkL/gLqCfMLIeZ3Y6/b/G5WGCbJ3vgP1B2w99tMPHxOHC3Zb6R6UREwN8NuXngzPWXwK1mdkPgG/5I8w9aEBX4smwB8KqZ5TSzq/F3A0uUmv00+L8c65PkgvmCZtY+MG8iUMvMbg/sa58ALjnbG5nZi2Z2eSBnJP4zIXuA1fiPTSfwXwebw8xewn+sS417zKymmeXBfw3QKHfKqKnOP9DQZKC/mRUIHLcrmdl1F7EOSZcrgH/01xecc6feiuFc2zQl+uIfsCAS/3Y7iv9viDz4j9FJbScwiNJZDANeMLPigWu2XsL/+yYhSIWPeKkT/m+7duC/Zqe5c+5oYN6lwA/4+2ovx79T6pj0xYEdTA/gRYBAIfEY/lF0Psb/h/tpnHO7gcb4v7n61cz24+8LvpfkFzuezSPAPwOvewn/afpzGYK/C95MYAP+b5oeD2TZG3i/T4HN+L9ZjDvz25y2HsfwX9jZBX8/7juBMed5WWP83RmSPtoFfg51zm1LfAD/xd+X+8YUxPnAko++szAl6yAikl6cczuBocCLgbMFrfFfW7IT/xmDf/C/v4PuAq7E32Xq5cDrEt/novfTgdePBf4FDA90pVqOf+AFnHPxQHv8f4jvwn+B/exzvR3wGf6zLlvwD/Bzi3PuAP5rlr7HX6j9hf9Yc9Zucyn0Bf4vJLfhHzX0ibMsdy+QE1iB/3g0irN3KT/XOiRVH/81sQOSHl/g3Ns0hSYGcj6I/9/6L/z/tivwXz+c1H+BmoFudf93hvd6HX/h/BuwDP/gCLqReIiy5JcHiIiIiGRv5r9NQWXn3D1eZ/GKmU3HP2jD2a6XFcl0dMZHRERERESyPBU+IiIiIiKS5amrm4iIiIiIZHk64yMiIiIiIlmeCh8REREREcnyMs39OYoVK+YqVKjgdYzTrF27lpMnT1KtWjWvo4iIpLuFCxfGO+dSepPdbCVUj1MiItnJuY5TmabwqVChAgsWLPA6xmn+/vtvChUqhP9myyIiWZuZ/eV1hlAVqscpEZHs5FzHqUxT+ISqwoULex1BRERERETOQ9f4pIH33nuPJ544282MRURERETEayp80kBcXBzr1q1DQ4OLiIiIiIQmdXVLA2+//bau8RERERGRTOP48ePExcVx5MgRr6NclMjISKKiooiIiEjxa1T4pIHEosc5pwJIREREREJeXFwc+fPnp0KFCpnu71fnHLt27SIuLo6KFSum+HXq6pZG3nrrLerXr6/ubiIiIiIS8o4cOULRokUzXdED/pMORYsWveCzVSp80kiFChW48sorOXr0qNdRRERERETOKzMWPYkuJrsKnzTSsWNHPv74YyIjI72OIiIiKWRmZc1smpmtNLPfzaxHYPorZrbZzJYEHjd7nVVEJKu5//77KVGiBLVr1w5O2717N82bN6dKlSo0b96cv//+O80+T4VPGtu5c6fXEUREJOVOAE8752oADYFHzaxmYN47zrm6gcd33kUUEcmaunTpwg8//JBsWt++fWnatClr166ladOm9O3bN80+T4VPGnrrrbcoW7YsBw8e9DqKiIikgHNuq3NuUeD5fmAlUCajc4wbN47vv/8+oz9WRMRT1157LUWKFEk27dtvv6Vz584AdO7cmf/7v/9Ls8/TqG5pqEWLFkRGRpKQkOB1FBERuUBmVgGoB/wKXAU8Zmb3AgvwnxU6rb+FmXUDugGUK1fuoj/79ddfJ3/+/Nx0000X/R4iIlnB9u3bKVWqFAClSpVix44dafbeOuOThho0aMBTTz1FgQIFvI4iIiIXwMzyAaOBJ51z+4CBQCWgLrAV6H+m1znnBjnnYpxzMcWLF7/oz4+OjmbZsmUX/XoRkdTy+XzExsYC/nv8+Hw+vvzySwAOHTqEz+djxIgRAOzduxefz8eYMWMAiI+Px+fzMX78eAC2bduW8SuQAip80tiRI0eYNm2ahrUWEckkzCwCf9HzlXNuDIBzbrtzLsE5dxIYDFyRnhmio6PZuXMn27dvT8+PEREJeSVLlmTr1q0AbN26lRIlSqTZe6urWxqLjY2le/furFy5kurVq3sdR0REzsH846H+F1jpnBuQZHop59zWQPM2YHl65oiOjgZg2bJllCxZMj0/SkTkjKZPnx58HhERkaydJ0+eZO2CBQsmaxcrVixZ+5JLLrnoHK1ateLzzz+nd+/efP7557Ru3fqi3+tUKnzS2G233Ua5cuWoUKGC11FEROT8rgI6AcvMbElg2nNARzOrCzjgT+Ch9AyRtPBp1qxZen6UiEjI6NixI9OnTyc+Pp6oqCheffVVevfuzR133MF///tfypUrx8iRI9Ps81T4pLGSJUty88263YOISGbgnJsFnOkueBk6fHWJEiUoWbKkrvMRkWxl2LBhZ5w+derUdPk8XeOTDrZt28Z//vMfDhw44HUUERHJJDTAgYhI+kpx4WNmQ8xsh5ktTzLtrHe2NrM+ZrbOzFab2Q1Jpjcws2WBee8H+ldnKStWrOCJJ55g9uzZXkcREZFMIjo6muXLl+uWCCIi6eRCzvjEAjeeYfppd7YO3PW6A1Ar8JqPzCw8sPxA/Pc8qBJ4nOk9M7Wrr76atWvXcsMNN5x/YREREfyFz5EjR/jjjz+8jiIikiWluPBxzs0Edqdw8dbAcOfcUefcBmAdcIWZlQIKOOd+cf7xnocCbS4wc8jLmTMnlStX9jqGiIhkIkkHOBARkbSXFtf4PGZmvwW6whUOTCsDbEqyTFxgWpnA81OnZzn79++ne/fujB071usoIiKSCdSsWRMzU+EjIpJOUlv4nO3O1me6bsedY/oZmVk3M1tgZgt27tyZyqgZK2/evMyePZs1a9Z4HUVERDKBPHnyULlyZRU+IiLpJFXDWTvngreYNrPBwIRAMw4om2TRKGBLYHrUGaaf7f0HAYMAYmJizloghaKwsDAWLVpEjhwaMVxERFJGI7uJSHaUkJBATEwMZcqUYcKECed/wUVK1RmfwDU7iZLe2Xoc0MHMcplZRfyDGMwL3AV7v5k1DIzmdi/wbWoyhLLEomfbtm0eJxERkcwgOjqadevWcejQIa+jiIhkmPfee48aNWqk++dcyHDWw4BfgGpmFmdmXYG3A0NT/wZcDzwF4Jz7HfgGWAH8ADzqnEscn7M78Cn+AQ/+AL5Pq5UJRcOGDSMqKorVq1d7HUVEREJcdHQ0zjlWrFjhdRQRkQwRFxfHxIkTeeCBB9L9s1LcD8s51/EMk/97juXfAN44w/QFQO2Ufm5m17RpU5599lmKFCnidRQREQlxSUd2i4mJ8TiNiEj6e/LJJ3n77bfZv39/un+WLkBJZyVKlOCNN06r/0RERE5TqVIlcufOret8RCTD+Xy+8y7TsmVLevXqFVy+S5cudOnShfj4eNq1a5ds2enTp5/3/SZMmECJEiVo0KBBipZPrbQYzlpSYO7cuYwaNcrrGCIiEsLCw8OpVasWv/32m9dRRETS3ezZsxk3bhwVKlSgQ4cO/PTTT9xzzz3p9nnmv49o6IuJiXELFizwOsZFa9myJWvXrmXVqlX4x3UQEcl8zGyhc059sM4grY5T999/PxMnTmT79u3nX1hE5CKtXLkyQwYUSKnp06fTr1+/CxrV7UzrcK7jlM74ZJCPPvqIhQsXqugREZFzio6OZseOHezYscPrKCIiWYoKnwxSrlw58uXLB0BmOcsmIiIZL+kAByIi2YXP50vXe/iACp8MtXnzZq655homTpzodRQREQlRKnxERNKHCp8MVKJECcLCwjh27JjXUUREJESVLFmS4sWLq/AREUljGs46A0VERDBjxgyvY4iISIiLjo5W4SMiksZ0xscDzjlGjx5NQkKC11FERCQERUdH8/vvv3Py5Emvo4hIFpaZrzu/mOwqfDwwZcoU2rVrx4gRI7yOIiIiISg6OppDhw6xfv16r6OISBYVGRnJrl27MmXx45xj165dREZGXtDr1NXNA82bN2f8+PHcfPPNXkcREZEQVKdOHQB+++03Kleu7HEaEcmKoqKiiIuLY+fOnV5HuSiRkZFERUVd0GtU+HjAzGjZsiUA+/fvJ2/evISF6eSbiIj41apVCzNj2bJl3H777V7HEZEsKCIigooVK3odI0Ppr20Pbdy4kRo1avDpp596HUVEREJInjx5qFSpkgY4EBFJQyp8PFS2bFlat25NvXr1vI4iIiIhRiO7iYikLXV185CZ8eGHH3odQ0REQlB0dDTffvsthw8fJnfu3F7HERHJ9HTGJwQcP36cl19+meHDh3sdRUREQkR0dDQnT55kxYoVXkcREckSVPiEgLCwMCZPnswvv/zidRQREQkR0dHRAOruJiKSRtTVLQSEh4czdepU8uTJ43UUEREJEZUrVyYyMlKFj4hIGtEZnxCRWPT89ddfDBw40OM0IiLitfDwcGrWrKnCR0QkjajwCTHvv/8+zz33HNu3b/c6ioiIeCw6OprffvvN6xgiIlmCCp8Q07dvX+bPn0/JkiW9jiIiIh6rU6cO27dvz7R3VhcRCSUqfEJMREQElStXBmDo0KHMnz/f40QiIuKVxAEOdNZHRCT1VPiEqMOHD/Pqq6/yzjvveB1FREQ80qBBAwDmzZvncRIRkcxPo7qFqNy5czN9+nR1eRMRycaKFClCtWrVdLsDEZE0oDM+Iaxs2bLkzJmTgwcP8uyzz3Lw4EGvI4mISAZr2LAhc+fOxTnndRQRkUxNhU8m8OuvvzJgwADmzJnjdRQREclgjRo1YufOnaxfv97rKCIimZoKn0ygSZMmrFq1iubNmwNw9OhRjxOJiEhGadiwIYC6u4mIpJIKn0yiUqVKACxatIhKlSrpACgikk3Url2bfPnyMXfuXK+jiIhkaip8MpmiRYtSp04dKlSo4HUUERHJAOHh4VxxxRX6wktEJJVU+GQy5cuX57vvvqNUqVI45xg9erQueBURyeIaNmzI0qVLOXTokNdRREQyLRU+mdjEiRNp164do0eP9jqKiIiko0aNGpGQkMCCBQu8jiIikmmp8MnEbrnlFsaNG0fbtm0B2Lt3r8eJREQkPWiAAxGR1FPhk4mZGbfeeitmxu7du6lVqxb//ve/vY4lIpJpmFlZM5tmZivN7Hcz6xGYXsTMppjZ2sDPwl7mLFasGJUrV9YAByIiqaDCJ4vInTs3d955J02bNvU6iohIZnICeNo5VwNoCDxqZjWB3sBU51wVYGqg7alGjRrxyy+/6LpOEZGLpMIni8idOzf9+/enfv36ALz88sv06tWLkydPepxMRCR0Oee2OucWBZ7vB1YCZYDWwOeBxT4H2ngSMImGDRuyfft2/vzzT6+jiIhkSip8siDnHLt372b37t2EhemfWEQkJcysAlAP+BUo6ZzbCv7iCCjhYTTAf8YHUHc3EZGLpL+KsyAz4z//+Q+DBw8GYP369TRq1Ihly5Z5nExEJDSZWT5gNPCkc27fBbyum5ktMLMFO3fuTL+AQHR0NHny5NEAByIiF0mFTxYWHh4OwObNm9m3bx9FihQBUPc3EZEkzCwCf9HzlXNuTGDydjMrFZhfCthxptc65wY552KcczHFixdP15w5cuTg8ssv1xkfEZGLpMInG7jmmmtYvnw5ZcqUAeCuu+6iZ8+eHqcSEfGemRnwX2Clc25AklnjgM6B552BbzM625k0atSIxYsXc/jwYa+jiIhkOip8sgn/sd1/ticqKorChf83MuvixYs1SpCIZFdXAZ2AJma2JPC4GegLNDeztUDzQNtzjRo14sSJEyxcuNDrKCIimU4OrwNIxgoLC6Nfv37B9ty5c2nUqBFffvkld999t4fJREQynnNuFmBnmR1y9wdIvJHp3Llzufrqqz1OIyKSueiMTzZXp04dBg4cSOvWrQH48ccf6devH0ePHvU4mYiInKpEiRJceumlGuBAROQiqPDJ5vLkycPDDz9Mvnz5ABg3bhwffPABERERAPzxxx8cP37cy4giIpJEw4YNdSNTEZGLoMJHknn//fdZsmQJYWFhOOe4+eabadeundexREQkoFGjRmzdupVNmzZ5HUVEJFNJceFjZkPMbIeZLU8yrYiZTTGztYGfhZPM62Nm68xstZndkGR6AzNbFpj3viVedS8ho1ChQoD/Rqj9+/enR48eABw6dIiYmBgmTJjgYToRkewt8Uam6u4mInJhLuSMTyxw4ynTegNTnXNVgKmBNmZWE+gA1Aq85iMzCw+8ZiDQDagSeJz6nhIiwsLCaNmyJU2aNAFg+/bt5M+fn7x58wL+G6M+++yzbN682cuYIiLZSp06dcidO7fu5yMicoFSXPg452YCu0+Z3Br4PPD8c6BNkunDnXNHnXMbgHXAFYGbwBVwzv3i/J2ThyZ5jYS4ihUrMm3aNK6//noA5s+fzzvvvMOJEycAWLhwIcOHD+fYsWNexhQRydIiIiKIiYnRGR8RkQuU2mt8SjrntgIEfpYITC8DJO18HBeYVibw/NTpkgndeeedxMfHU758eQCGDh1Kt27dgvcMmjVrlg7MIiLpoFGjRixatIgjR454HUVEJNNIr8ENznTdjjvH9DO/iVk3M1tgZgt27tyZZuEk7RQoUCD4fMCAASxYsCA4ItxLL73EI488Epw/ZswYdc0QEUkDDRs25Pjx4yxevNjrKCIimUZqC5/tge5rBH7uCEyPA8omWS4K2BKYHnWG6WfknBvknItxzsUUL148lVElvYWHh1O1atVge9SoUcTGxgbbTz31FO+++26w/c477zBnzpwMTCgikjUk3shUZ9VFRFIutYXPOKBz4Hln4Nsk0zuYWS4zq4h/EIN5ge5w+82sYWA0t3uTvEaymCJFinDZZZcF20uWLOFf//oXAEeOHOG5555j0qRJAJw4cYKOHTsydepUT7KKiGQmpUqVonz58ip8REQuwIUMZz0M+AWoZmZxZtYV6As0N7O1QPNAG+fc78A3wArgB+BR51xC4K26A5/iH/DgD+D7NFoXCXGFCxcOXg8UGRnJ33//Tc+ePQHYunUr8+bNY/v27QD89ddfVKpUKVgYHT58mC1btuiGfSIiAY0bN2b27NnaL4qIpFCOlC7onOt4lllNz7L8G8AbZ5i+AKid0s+VrCsyMpLIyEgAypYtyx9//BE8gB87doz69etzySWXADB79myaN2/OtGnT8Pl8rF27lilTptChQweKFCni2TqIiHilSZMmDBs2jBUrVlCrVi2v44iIhLz0GtxA5KIkjghXpUoVRo4cGewqV61aNd5//33q1KkDwPTp03n00Uc5cOAAAN988w3XXXcd8fHxAGzevJm1a9dy8uRJD9ZCRCT9NW/eHIDJkyd7nEREJHNQ4SOZQtmyZXn88ceDZ3ceeOABNm3aRNmy/jE0EgumQoUKAfDxxx9To0YNjh8/DsDXX39Njx49gmeUdu3apWFgRSRTK1++PNWqVWPKlCleRxERyRRU+EimZGZERUUFC5727dszY8YMcuTw99685557GDZsGLly5QJgxYoV/PTTT8Hln3766WQj0H3yySe8/fbbwfa2bds4ePBgRq2OiMhFadGiBdOnT+fo0aNeRxERCXkqfCRLqlatGu3btw+2X3/9dZYtWxZs33333bz66qvB9syZM/n++++TzU/sRgLw/PPP069fv2B74cKFbNiwIb3ii4ikSIsWLTh8+DCzZ8/2OoqISMhL8eAGIllJ0qIG4Kuvvko2MtJTTz1FQkJCsL1ixQqS3kuqY8eO1KtXjxEjRgD+Pz6uvvpqXnrpJQA++ugjateuzbXXXgvAjh07KFKkSPCMlIhIWvD5fERERDB58mSaNGnidRwRkZCmMz4iAYnd4ABatmxJ69atg+2xY8cyaNCgYDs2NpbevXsH21FRURQtWhQA5xy9evXi22+/DbbLlStHnz59gu3mzZszbNgwAE6ePMnAgQNZvnx5sL1jxw5OnDiRTmsqIllFvnz5aNy4sQY4EBFJARU+IhehcePG1KtXL9geMmQIjz76KOAvoHbu3MnLL78M+AuZAQMG0KZNG8B/T6IjR45w7NgxAP7++28eeeSR4M1b4+PjKVmyJJ988gkA27dv58orr+S7774DYPfu3bz++uusXLkSgEOHDjFv3jz27t2b/isuIiGnRYsWLF68mB07dngdRUQkpKnwEUkHefPmpUCBAgCEh4fzyCOPcNVVVwGQJ08efv75Zzp37gz4b+y6efPmYDsyMpL3338/2E3u2LFjFCpUKDhQw6ZNm3jxxRdZtWoV4O+Gd+WVVzJz5kwAfv31VwoWLMhPP/0EwG+//cZtt93GihUrANiwYQPvvvsu27ZtA/wj3C1YsIDDhw8D/jNSuiGiSObRokULAH788UePk4iIhDYVPiIeCwsLo3Tp0sGhuAsUKMDjjz9OdHQ04B/Ke9KkSTRt6r9X8GWXXcbRo0dp2bIlAJUqVWL8+PFcccUVABQtWpT77rsvONT3/v37WbduXbDr3JIlS3jqqaeChc/UqVO5/PLLWb9+PUBwNLx169YBMH78eK699lq2b98OwKxZs+jVqxf79+8H/IXX8OHDg2ewduzYwZo1a4L3UFIRJZK+6tWrR5EiRTSstYjIeajwEcmEcubMSUREBOA/Y9SyZUtKliwJQOXKlXn33XepUqUKAFdddRXLli0L3vy1VatW7N69m9q1awNwzTXX8O2331K+fHkAqlevTs+ePSlWrBjg77oXHh5Ozpw5AVi2bBkDBw4MFjbjx4+nY8eOwcJq0KBBVKtWLTg4xMsvv0y+fPmCBdB7770XPJsF/oElnnjiiWB78uTJwW5+AEuXLmXGjBnB9vbt29m6dWvqN6JIFhEeHk6zZs2YPHmyvmgQETkHFT4i2Ux4eDiFCxcOjjBXqlQpWrVqRb58+QCoX78+ffv2DZ6BatmyJdOmTaNw4cIAdO/enYMHD1KwYEEAHnroIVasWEHu3LkBuO222/jqq6+ChdlVV13FE088ERw8Im/evMGiCmDVqlVMnz492B4xYgSvvfZasN2vXz/uu+++YLtHjx74fL5gu127dlx++eXB9sMPP8xdd90VbL/00ku88MILwfaHH36YbKCKUaNGJRvKfNasWSxZsiTY/uOPP4Jnx8B/jVbSEf9EQkGLFi3YsmVLsEuriIicToWPiKRKoUKFqFGjRrCwqVWrVrLC44YbbuDNN98Mth944AHGjBkTbL/22mv89ttvwfbAgQP5/fffg+1XX32VUaNGBdvdu3fnrbfeCrZbt25Np06dgu2oqCjKlSsXbG/ZsoXNmzcH22PHjmXcuHHB9ptvvsnAgQOD7YcffpjXX3892L7pppvo2bNnsF2jRg26du0abEdHR/PMM88E202bNk12z6cOHTrw+eefB9tPPvkkEydOBPzdAPv27cucOXMASEhI4Isvvgj+8XrixAmmTp1KXFxcsL1y5Ur27NkD+AfO+Pvvvzl+/DiSvSUO0a/R3UREziHxQuZQfzRo0MCJiKS1Xbt2ufj4+GB76dKlbuXKlcH2+PHj3axZs4LtDz74wI0fPz7YfuaZZ9xXX30VbLdv3959/PHHwXZMTIzr37+/c865kydPuqJFi7rXX3/dOefcsWPHHOBee+0155xzBw4ccIB7++23nXPO7d692wHu3Xffdc45t23bNge4Dz/80Dnn3MaNGx3gPv30U+ecc+vWrXO5cuVyX3/9tXPOuTVr1rhy5cq5CRMmpHYzBQELXAgcE0Lx4fVxqnr16u7GG2/0NIOIiNfOdZzS3RRFJFsrUqRIsnbitVCJEgeRSJQ4bHmif/3rX8na33zzTbL2/Pnzg8/NjPj4+GA7R44cHD58OHi2LHfu3KxduzaYKV++fMyYMYNLL70U8A98MXz4cOrXrw9AwYIFeffdd2nYsGFw/pNPPkn16tWD73f99ddTokSJ820GyQJatGjB4MGDOXLkCJGRkV7HEREJOeYyyYWQMTExbsGCBV7HEBHJ1sxsoXMuxuscocjr49TEiRNp2bIlP/74Y3AUSBGR7OZcxyld4yMiIpIFXHfddURERGhYaxGRs1DhIyIikgXky5ePq666SgMciIichQofERGRLKJ58+YsXryYHTt2eB1FRCTkqPARERHJIlq0aAHAjz/+6HESEZHQk21GdUt6w8OzadmyJb169Qou36VLF7p06UJ8fDzt2rU77+tPXf7pp5/m1ltvZfXq1Tz00EPnff2py7/55ps0btyYOXPm8Nxzz5339acu/8knn1CtWjXGjx9P//79z/v6U5cfNWoUxYoVIzY2ltjY2PO+/tTlE29K2a9fPyZMmHDe1ydd/pdffmH06NEA9OnTh19++eWcry1atGiy5Xft2hW8SWW3bt1Ys2bNOV9ftWrVZMsXLVo0eK+Ytm3bsmvXrnO+vlGjRsmWb9SoUbLfpfPR755+9xKXT+/fvaQ3i5Wsp169ehQtWpTJkycnu5+WiIjojI+IiEiWER4eTrNmzZg8eTKZZdRWEZGMouGsRUQkxTSc9dmFynFqyJAhdO3alWXLllG7dm2v44iIZCgNZy0iIpJNNG/eHIAffvjB4yQiIqFFhY+IiEgWUrZsWerXr8+oUaO8jiIiElJU+IiIiGQx7du359dff2Xjxo1eRxERCRkqfERERLKY9u3bA+isj4hIEip8REQyof3797N79+5g+6+//mLdunXB9sKFC0l6of2kSZOYNWtWhmYU71SqVIl69eoxcuRIr6OIiIQMFT4iIilw8OBBNm/eHBwiePPmzcyePTs4f/ny5YwYMSLYnjFjBu+++26wPWbMGPr06RNsDxo0iK5duwbbr7/+Om3atAm2e/TowXXXXRdsd+jQgcsuuyzY7tixY/AidoAHH3yQTp06BdtPP/108H4+AK+++ioffvjhha62ZGLt27dn7ty56u4mIhKgwkdEMg3nXLDwOHz4MGvXruXw4cMA7NixgwkTJrB3714A1q5dyzvvvBO8Aej8+fN59NFH2blzJwCTJ0/mhhtuCLaHDx9O1apViY+PB+CDDz4gR44c7NmzJ9iOioriyJEjgH/I4KuvvpoTJ04AMGLECDp27BjM991339G7d+9g9nnz5vH1118H29u2bUt2hiZ37tzky5cv2K5VqxaNGjUKtlu3bp2sUHr44Yd55plngu0XX3yRN998M9h+//33+c9//hNsDx8+PFkhJlmfuruJiJwi8Q+JUH80aNDAiUjmcfz4cbdhwwa3Z88e55xzBw4ccGPHjnV//vmnc865HTt2uFdffdUtW7bMOefc+vXr3R133OF+/fVX55xzS5cudbVq1XIzZ850zjk3Y8YMFxYW5qZNm+acc27y5MkOcLNmzXLOOTd+/HgHuHnz5jnnnBszZowD3OLFi51zzo0ePdoVLVrUrVmzxjnn3IQJE9yVV17p4uLinHPOTZkyxXXs2NHt2rXLOefc7Nmz3XPPPecOHDjgnHNu8eLFbtCgQe7o0aPOOefWrl3rJk2a5E6cOOGcc27btm1u5cqV7uTJk8455w4dOhR8bVYCLHAhcEwIxUcoHqfq1avnGjZs6HUMEZEMc67jlM74iEjQtm3bgteNOOcYO3Ysy5YtA+DEiRP07t2bKVOmAHDo0CGaNWvG8OHDAdi1axeFChVi4MCBgP8MTMWKFYPzd+7cyW233cZPP/0EwJ49e3j55ZdZunQpAMePH2fp0qXBMyx58+alevXq5MmTB4Dy5cvz3HPPERUVBUB0dDRffPEFVapUAeCqq67i119/pWbNmgDcfPPN7Nmzhzp16gBw++23Ex8fH1z+lltuYe7cuZQpUwaAZs2a8fXXX1OkSBEAGjduzBtvvEHevHkBqFu3Lg8++CA5c+YEoHLlyrRo0YLw8HAASpYsSfXq1TEzwH8GJ/G1Il5J7O62adMmr6OIiHhOhY9IFrJ161Y2b94cbI8cOZKpU6cG2z169GDQoEHBdr169Xj++eeD7UqVKvHGG28AYGZ06NCBr776CoDw8HDef/995s+fD0DOnDk5fPgwCQkJAOTLl4/OnTsHC48iRYowZMgQrr/+egBKly7NokWLgtexVKpUiePHj3P33XcDULVqVVatWkWLFi2C80eNGkWDBg0Af+Hz2muvUblyZQAuueQS7rnnHkqUKAFA4cKFueKKK4LFRq5cuShYsCBhYdrNydmZ2RAz22Fmy5NMe8XMNpvZksDjZi8zpoa6u4mI/E8OrwOISHIJCQnBswiLFi3i0KFDXH311QC89957JCQk0LNnTwDatWtHvnz5iI2NBaBFixZUqVKFMWPGAP7rPqKjo2natCkAixcvJjIyMvhZ1157LdWrVw+2P/roo2TtBQsWULJkScBfCB06dCg4L0eOHMku7s+VKxfvvfdesB0ZGcl9990XbOfMmZN69eoF22FhYSpKJBTEAh8AQ0+Z/o5zrl/Gx0lblStXpm7duowcOZKnnnrK6zgiIp5S4SOSzk6ePBn8A3/p0qVs2rSJli1bAvDhhx+yevVq3n//fcA/cteaNWtYtGgRAC+88AI7duwIDks8ffr0ZIVPdHQ0uXPnDn7WG2+8QcGCBYPtqVOnJrtgfubMmcmyJS1UADp37pysHR0dffErLpIJOOdmmlkFr3Okp/bt2/P888+zadMmypYt63UcERHP6OtWkVTYuXMnc+fODbbHjRvHI488Emw/88wzXHLJJcH2wIEDuf/++4PtjRs3snLlymC7VatWdOnSJdju378/X375ZbA9duxYxo0bF2y//PLLyUb2atWqVbIhkMuUKZOsEBKRFHvMzH4LdIUrfLaFzKybmS0wswWJIwSGGnV3ExHxU+Ejcg4bNmzg888/Dw5hPGLECOrWrRvs8jV48GAaNWoUHFJ55cqVjB8/nuPHjwP+rmTdu3cPDnHcp08ffv755+D7/+tf/woOFgBw11138cQTTwTbNWrUSNb1TEQyxECgElAX2Ar0P9uCzrlBzrkY51xM8eLFMyjehalSpUqwu5uISHamwkeytU2bNvHhhx+yY8cOACZMmEDJkiVZs2YNAD///DNdunQJ3gAwX758lCtXjoMHDwL+b1K///774DU5zz77LJs2bSIiIgKAli1b8uqrrwZH+ipfvjzVqlXL0HUUkQvjnNvunEtwzp0EBgNXeJ0ptdq3b88vv/yi0d1EJFtT4SNZ2p49exg2bBh//fUX4L9Yv0yZMkyfPh2AdevW8dhjj/Hbb78BUK5cOdq0aRMcsvjWW29l7dq1XHrppYB/CORx48aR+M1ulSpVuPHGG4PLi0jmZ2alkjRvA5afbdnMIrG72+jRoz1OIiLiHRU+kqk559iwYQNbt24FYPv27fh8PsaOHQv47y1z1113Be8dU7p0aVq0aBG87qVRo0Zs2bKFJk2aAFCnTh0++eQTKlSoAPiHSK5cuTI5cmgcEJGsyMyGAb8A1cwszsy6Am+b2TIz+w24Hsj0w6FVqVKFyy67TN3dRCRbU+EjmYpzjv79+/P9998D/ptoVqpUicGDBwP+QuXEiRPBa2oqVKjAsmXL6NChA+AvfD777LPgsMqRkZGUKlVKwyqLZFPOuY7OuVLOuQjnXJRz7r/OuU7OuWjnXB3nXCvn3Favc6aF9u3bM2fOHHV3E5FsS3/tSchJSEhg165dwfadd95Jnz59AP+9ZAYMGMD48eMByJs3L1988QV33HEH4L9XzKxZs7j99tsB/003a9eunWzIZxGR7Ejd3UQku1P/HfHc5s2b2bhxI40aNQL8I6Hlzp2bH3/8EYBChQpRoECB4PKrVq0if/78wfbdd9+dsYFFRDKhqlWrUqdOHUaOHMmTTz7pdRwRkQynwkcy3PLly/n111/p2rUrAL169WL27NnBkdN69OiRbPlPPvkkWTtp0SMiIil3xx138MILL7BhwwYqVqzodRwRkQylrm6S7lavXs2bb75JQkICAN988w0PPfRQ8F44vXv3ZvTo0cHrcu64445g1zURkTMxs5fO8XjR63yhqnPnzoSFhfHpp596HUVEJMOlSeFjZn8GRsBZYmYLAtOKmNkUM1sb+Fk4yfJ9zGydma02sxvSIoOEjq1bt9KvXz8S72K+cOFCXnjhBVavXg3AY489xpYtW8iTJw8Al112GZdffnnwXjciIilw8AwPB3QFnvUwV0iLiorilltuYciQIcEbLYuIZBdpecbneudcXedcTKDdG5jqnKsCTA20MbOaQAegFnAj8JGZhadhDslghw8fZuTIkcHCZsuWLfzjH/9g9uzZALRu3Zpdu3ZRs2ZNAEqUKEGJEiU8yysimZ9zrn/iAxgE5AbuB4YDl3oaLsR169aNbdu2BQeJERHJLtKzq1tr4PPA88+BNkmmD3fOHXXObQDWkQXuip3d7Nixg3Xr1gFw9OhROnbsyLBhwwCoV68eGzdupE2bNoB/5LXChQuf7a1ERC5KoGfB68Bv+K9Zre+ce9Y5t8PjaCHtxhtvJCoq6rTrJ0VEsrq0KnwcMNnMFppZt8C0kon3Pgj8TPyKvwyQ9CYCcYFppzGzbma2wMwWJHabEu8cOXIE8N9Lp0GDBvTu3Rvwj7q2aNEiXnjhBQDCwsIoW7asZzlFJOszs38D84H9QLRz7hXn3N8ex8oUcuTIwQMPPMDkyZPZsGGD13FERDJMWhU+Vznn6gM3AY+a2bXnWPZMF3K4My3onBvknItxzsUUL148LXLKRXriiSdo2LAh4L+XzsCBA3n55ZeD8+vUqUOOHBokUEQyzNNAaeAFYIuZ7Qs89pvZPo+zhbyuXbtqkAMRyXbSpPBxzm0J/NwBjMXfdW27mZUCCPxM7HoQByQ9HRAFbEmLHJJ2Vq5cSY8ePTh27BgAV111FbfffntwZLaWLVsSHR3tZUQRycacc2HOudzOufzOuQJJHvmdcwXO/w7ZmwY5EJHsKNWFj5nlNbP8ic+BFsByYBzQObBYZ+DbwPNxQAczy2VmFYEqwLzU5pDU27dvHwcOHADgzz//ZPDgwSxduhSAO++8k5deeonwcI1DISLeM7NGpqEgU0WDHIhIdpMWZ3xKArPMbCn+Amaic+4HoC/Q3MzWAs0DbZxzvwPfACuAH4BHnXMJaZBDUmH79u2ULVuWgQMHAtCiRQu2bdvG5Zdf7nEyEZEz6gwsMrPhZtbFzC7xOlBmo0EORCS7SfVFGc659cBlZ5i+C2h6lte8AbyR2s+W1Fm6dCkrV66kQ4cOlCxZkmeeeYZmzZoBEB4eToEC6i0iIqHJOfcwgJlVx399aayZFQSm4f9Sbba+VDu3xEEOXnnlFTZs2EDFihW9jiQikq7SczhrCXH/+te/ePrpp4P9u59//nnq1avncSoRkZRzzq1yzr2D/1YJTYBZQHvgV0+DZRKJgxwMHjzY6ygiIulOhU828tdff9GhQwfi4uIAePvtt1m+fDkREREeJxMRSbXRQHPn3HfOuceT3ExbzkGDHIhIdqLCJxtxzjF16lSWLFkC+A94urGoiGQRrYFyZvaVmVX1Okxm0q1bN7Zv365BDkQky1Phk8X169ePxx57DIAKFSqwadMmWrZs6XEqEZG05ZxLcM59ADwGdDMzXUeaQhrkQESyCxU+WdzOnTvZunVr8P47kZGRHicSEUl7ZtbSzPoA7+K/P1wZbxNlHomDHEyePJkNGzZ4HUdEJN2o8MliDhw4wKOPPsrixYsBeOuttxg9erTuvyMiWV0hYCLwgHOug3Oui7dxMpfEQQ4GDRrkdRQRkXSjwieLOXbsGOPGjWPmzJkAhIXpn1hEsoXbgQeAu8ws2sy087sAUVFRtGnThoEDB7Jv3z6v44iIpAsdGLKAw4cP8/HHH+Oco0iRIqxYsYIePXp4HUtEJN2YWbmkbefc7cDbwB78w1l/40GsTK1Pnz7s3bs3eCNrEZGsRoVPFvDNN9/QvXt3Zs2aBUD+/Pk9TiQiku5+MLN4M/vZzD4ys4eB8sA059xLzrl2XgfMbGJiYmjRogUDBgzg8OHDXscREUlzKnwysUOHDgHQqVMnfv31V6655hqPE4mIZAznXE2gNPAEMBeoDLwIrDIzXaF/kZ577jl27NjBZ5995nUUEZE0p8Ink/rss8+oUaMGW7ZsISwsjCuuuMLrSCIiGco5d8w5txgYC/wKbAMOA0s9DZaJXXvttTRu3Ji3335bNzQVkSxHhU8mVb9+fa655hoKFizodRQRkQxnZtXMrKeZ/QTMARoBXwE1nHNtPA2XiZkZffr04a+//mLYsGFexxERSVPmnPM6Q4rExMS4BQsWeB3DUwcOHGDSpEm0bdvW6ygikk2Z2ULnXEwI5DgJLAb6AuOcc0c9jpRljlPOOerWrcuxY8f4/fffNTqoiGQq5zpOaW+Wibz99tt06NCB9evXex1FRMRr3YFfgMeATWa20sy+MbMXzayNt9Eyt8SzPqtWreL//u//vI4jIpJmdMYnEzly5Ajz5s3j2muv9TqKiGRToXLG51RmFgXUAaKB2s65ThmdISsdpxISEqhevToFCxZk/vz5mJnXkUREUkRnfDKxPXv28Oijj7J//34iIyNV9IiInIFzLs45951z7l9eFD1ZTXh4OM8++ywLFy5kypQpXscREUkTKnxC3OzZs4mNjeW3337zOoqIiGQjnTp1okyZMrz55pteRxERSRMqfELcLbfcwoYNG7jqqqu8jiIiItlIrly56NWrFzNmzGDOnDlexxERSTUVPiHIOceLL77Izz//DECJEiU8TiQiItnRgw8+SNGiRXnrrbe8jiIikmoqfELQvn37GDlyJN9++63XUUREJBvLmzcvTz75JBMmTGDRokVexxERSRUVPiGoYMGCzJ07l7ffftvrKCIiks09/vjjFCtWjF69epFZRoIVETkTFT4h5Oeff+aZZ57h5MmTFCpUSDeNExERzxUsWJBXX32VadOmMX78eK/jiIhcNP1lHUK+//57vv32Ww4cOOB1FBERkaBu3bpRvXp1evXqxbFjx7yOIyJyUVT4hJA333yTBQsWUKBAAa+jiIiIBOXIkYN+/fqxdu1aPv74Y6/jiIhcFBU+IWDx4sVs2LABgPz583ucRkRE5HQ333wzzZo149VXX+Xvv//2Oo6IyAVT4eMx5xzdunWjVatWumhURERClpnRv39//v77b15//XWv44iIXDAVPh4zM0aPHs2QIUMwM6/jiIiInFWdOnW4//77+c9//sO6deu8jiMickFU+HgocRCDcuXKcfnll3ucRkRE5Pxee+01cubMSe/evb2OIiJyQVT4eOT48eNcffXVPP30015HERERSbFSpUrRu3dvRo8ezc8//+x1HBGRFFPh4xHnHG3atOG6667zOoqISLZlZkPMbIeZLU8yrYiZTTGztYGfhb3MGIp69uxJVFQUPXv25OTJk17HERFJERU+HsmZMyevvPIKrVq18jqKiEh2FgvceMq03sBU51wVYGqgLUnkyZMneAuGYcOGeR1HRCRFVPhksJMnT3Lfffcxf/58r6OIiGR7zrmZwO5TJrcGPg88/xxok5GZMou7776bmJgY/vGPf7Bnzx6v44iInJcKnwz2119/MWXKFJYvX37+hUVExAslnXNbAQI/S2TEh/p8PmJjYwH/daA+n48vv/wSgEOHDuHz+RgxYgQAe/fuxefzMWbMGADi4+Px+XyMHz8egG3btuHz+fjhhx8A2LRpEz6fjx9//BGA9evX4/P5mDFjBgCrV6/G5/MxZ84cAJYvX47P5wt+SbdkyRJ8Ph9LliwBYP78+TRp0oRevXqxY8cOOnXqhM/nY/Xq1QDMmDEDn8/H+vXrAfjxxx/x+Xxs2rQJgB9++AGfz8e2bdsAGD9+PD6fj/j4eADGjBmDz+dj7969AIwYMQKfz8ehQ4cA+PLLL/H5fBw/fhyA2NhYfD5fcFsOHjyYZs2aBdsfffQRN910U7D93nvvJetx0a9fP9q2bRts9+3blw4dOgTbr732Gvfcc0+w/dJLL3HfffcF23369KFbt27Bdq9evXj00UeD7SeffJInn3wy2H700Ufp1atXsN2tWzf69OkTbN9333289NJLwfY999zDa6+9Fmx36NCBvn37Bttt27alX79+wXarVq147733gu2bbrqJjz76KNhu1qwZgwcPDrYz2++ez+cL/h01Z84c/e5lsd+99JQjXd9dTlOxYkU2bNigoatFRLIAM+sGdAP/CJ3ZTa1atfjHP/5B3759qVOnjtdxRETOyTLLTTNjYmLcggULvI6RKvv27SN//vwqekQk0zKzhc65GK9zpCUzqwBMcM7VDrRXAz7n3FYzKwVMd85VO9/7ZIXj1MU4cuQIdevW5ciRIyxbtoz8+fN7HUlEsrFzHafU1S0DPfTQQ1x99dVklmJTRCSbGgd0DjzvDHzrYZaQFxkZyZAhQ9i4cWOyLjMiIqFGhU8GatmyJR06dNAZHxGREGFmw4BfgGpmFmdmXYG+QHMzWws0D7TlHBo3bkyPHj348MMPmTlzptdxRETOSF3dREQkxbJiV7e0kt2PUwcPHqROnTqEhYWxdOlS8uTJ43UkEcmG1NXNY3v37iU2NpajR496HUVERCRd5M2bl08//ZR169YlGxVKRCRUqPDJAN988w333XcfK1as8DqKiIhIurn++ut56KGHeOedd5g7d67XcUREklHhkwEeeOAB5s6dS7169byOIiIikq7efvttSpcuzf3336+eDiISUlT4ZAAz48orr/Q6hoiISLorUKAAgwcPZuXKlfTu3dvrOCIiQSp80pFzjpYtWzJ06FCvo4iIiGSYG2+8kccff5x3332X0aNHex1HRATwsPAxsxvNbLWZrTOzLPmV0J49ezhy5AjHjx/3OoqIiEiG+ve//83ll1/O/fffz7p167yOIyLiTeFjZuHAh8BNQE2go5nV9CJLeipcuDA//vgj999/v9dRREREMlSuXLkYOXIk4eHhtG/fnsOHD3sdSUSyOa/O+FwBrHPOrXfOHQOGA609ypIutm3bxt69ewF0w1IREcmWypcvzxdffMGSJUvo0aOH13FEJJvL4dHnlgE2JWnHAel69b/P50vPtz/N2rVr2blzJw0bNiQsTJdSiUhomD59utcRJJu55ZZb6N27N3379uWaa66hU6dOXkcSkWzKq7/Iz3QKxJ22kFk3M1tgZgt27tyZAbHSTsmSJbn00ktV9IiISLb32muvcd111/Hwww/z+++/ex1HRLIpr874xAFlk7SjgC2nLuScGwQMAoiJiTmtMLoQ+pZTRETEGzly5GDYsGHUq1ePdu3aMX/+fPLly+d1LBHJZrw6HTEfqGJmFc0sJ9ABGOdRljQ3f/58faMlIiKSRKlSpRg2bBhr1qzhwQcfxLlUfZ8pInLBPCl8nHMngMeAScBK4BvnXJapFHr37k2HDh28jiEiIhJSrr/+el5//XWGDx/OP//5T6/jiEg241VXN5xz3wHfefX56enzzz9ny5bTeu6JiIhke71792bNmjW88sorlC9fni5dungdSUSyCc8Kn6wsKiqKqKgor2OIiIiEHDNj0KBBbN68mQcffJAyZcrQvHlzr2OJSDagIcfS2BdffMHEiRO9jiEiIhKyIiIiGDVqFDVr1qRt27YsXbrU60gikg2o8Elj//73vxkyZIjXMUREREJagQIFmDhxIgUKFOCWW24hLi7O60giksWp8EljCxcuZODAgV7HEBERCXlRUVF8//337N+/n5tvvpm9e/d6HUlEsjAVPmksIiKCEiVKeB1DREQkU4iOjmbMmDGsXLmSdu3acezYMa8jiUgWpcInDT355JPExsZ6HUNERCRTadq0Kf/973/58ccfuffeezlx4oTXkUQkC9KobmkkISGBuXPnkidPHq+jiIiIZDr33nsvO3bs4B//+AdmxhdffEGOHPozRUTSjvYoaSQ8PJy5c+dy8uRJr6OIiIhkSr169QLgH//4B845vvzySxU/IpJmtDdJY2Fh6j0oIiJysXr16oWZ0atXL5xzfPXVVyp+RCRN6K/0NOCco0GDBnz00UdeRxEREcn0nn76afr3788333zDXXfdxfHjx72OJCJZgL5CSQP79++nWrVqFC1a1OsoIiIiWULPnj0xM3r27Ilzjq+//pqIiAivY4lIJqbCJw0UKFCAr7/+2usYIiIiWcpTTz0FkKz4yZkzp8epRCSzUle3NLBv3z6vI4iIiGRJTz31FO+++y6jR4/mpptu0k1OReSiqfBJpT179lC8eHE+/PBDr6OIiIhkST169GDo0KHMnDmTq6++mk2bNnkdSUQyIRU+qZSQkMBzzz3HVVdd5XUUERGRLKtTp0788MMPbNy4kYYNG7JkyRKvI4lIJqPCJ5WKFi3Kyy+/TN26db2OIiIikqU1bdqUWbNmERYWxjXXXMOkSZO8jiQimYgKn1RwzjFv3jwSEhK8jiIiIpItREdHM3fuXC699FJuueUWhgwZ4nUkEckkVPikwtq1a7nyyisZOnSo11FERESyjTJlyvDzzz/TpEkTunbtyvPPP68vIUXkvFT4pELJkiUZMWIEzZo18zqKiIhItlKgQAEmTpxI165defPNN7n55puJj4/3OpaIhDAVPqlQsGBB7rjjDsqWLet1FBERkWwnIiKCwYMHM2jQIKZPn079+vWZN2+e17FEJESp8EmFhQsXsnbtWq9jiIiIZFtmxoMPPsjs2bMJCwvj6quvZuDAgTjnvI4mIiFGhU8qdO/ene7du3sdQ0REJNuLiYlh0aJFNGvWjEceeYTOnTtz6NAhr2OJSAjJ4XWAzGzQoEEcP37c6xgiIiICFClShAkTJvDGG2/w8ssvs2TJEkaMGEGNGjW8jiYiIUBnfFKhbt26XH755V7HEBERkYCwsDBefPFFfvjhB7Zu3Uq9evV45513OHnypNfRRMRjKnwu0p9//snEiRN1Gl1ERCQEtWjRguXLl9OiRQt69uxJ06ZN+fPPP72OJSIeUuFzkcaPH0/Lli3Zu3ev11FERETkDEqWLMm3337LkCFDWLhwIXXq1OGzzz7TwAci2ZQKn4vUuXNn5syZwyWXXOJ1FBERETkLM+O+++7jt99+o379+tx///20adOG7du3ex1NRDKYCp+LVKBAARo1aoSZeR1FREREzqNChQr89NNPDBgwgEmTJlGzZk3++9//6tofkWxEhc9FGjRoEMuXL/c6hoiIpBMz+9PMlpnZEjNb4HUeSb2wsDCeeuopFi9eTK1atXjggQe45pprWLZsmdfRRCQDqPC5CHv27OGhhx7iu+++8zqKiIikr+udc3WdczFeB5G0U6NGDWbMmMFnn33GmjVrqFevHr169eLAgQNeRxORdKTC5yIULFiQrVu30rVrV6+jiIiIyEUwM7p06cKqVau4//776d+/PzVq1GDMmDEa/EAki1LhcxHMjEsuuYSiRYt6HUVERNKPAyab2UIz6+Z1GEkfRYsWZdCgQcyZM4ciRYrQtm1bWrRowZIlS7yOJiJpTIXPRRg/fjyfffaZ1zFERCR9XeWcqw/cBDxqZteeuoCZdTOzBWa2YOfOnRmfUNJMo0aNWLhwIe+99x6LFi2ifv36dO7cmY0bN3odTUTSiAqfixAbG8uAAQO8jiEiIunIObcl8HMHMBa44gzLDHLOxTjnYooXL57RESWN5ciRgyeeeII//viDZ555hhEjRlC1alV69+7Nnj17vI4nIqmkwucijBw5khkzZngdQ0RE0omZ5TWz/InPgRaAhvLMJgoVKkTfvn1Zs2YNd955J2+//TaVK1fm3Xff5ciRI17HE5GLpMLnIoSFhVGkSBGvY4iISPopCcwys6XAPGCic+4HjzNJBitXrhyff/45CxcupF69ejz11FNceumlvPfeexw+fNjreCJygVT4XKANGzbw7LPPsn79eq+jiIhIOnHOrXfOXRZ41HLOveF1JvFOvXr1mDJlCtOmTaN69eo8+eSTVKxYkf79+3Pw4EGv44lICqnwuUCrVq3i3XffZd++fV5HERERkQzk8/n46aefmDlzJtHR0fTq1YuKFSvy9ttv6x5AIpmACp8LdNNNN3Hw4EHq1KnjdRQRERHxwDXXXMOUKVOYPXs29evX59lnn6Vs2bL07t2buLg4r+OJyFmo8LkIOXLkICxMm05ERCQ7a9y4MT/88AO//vorzZs359///jcVK1bknnvuYeHChV7HE5FT6K/3C9SrVy++/PJLr2OIiIhIiLjiiiv45ptv+OOPP3j88ccZN24cMTEx+Hw+vv32WxISEryOKCKo8Lkgzjl+/PFHli1b5nUUERERCTEVKlRgwIABbNq0if79+7NhwwbatGnDpZdeymuvvcaWLVu8jiiSrZlzzusMKRITE+MWLFjgdQzAXwCZmdcxREQynJktdM7FeJ0jFKX2OOXz+c67TMuWLenVq1dw+S5dutClSxfi4+Np167deV9/6vJPP/00t956K6tXr+ahhx467+tPXf7NN9+kcePGzJkzh+eee+68rz91+U8++YRq1aoxfvx4+vfvf97Xn7r8qFGjKFasGLGxscTGxp739acuP336dAD69evHhAkTzvv6pMv/8ssvjB49GoA+ffrwyy+/nLa8c474+Hi2bNnCnj17CA8Pp1WrVkRERFCgQAEGDx4MQLdu3VizZs05P7tq1aoMGjQouHzRokV56623AGjbti27du065+sbNWqUbPlGjRol+106H/3uZa7fvaSKFi2abPldu3Yl+10Ktd+9xHW9WOc6TuVI1TtnUyp6RERE5HzMjOLFi1O8eHFy5crFZZddxmeffUZ8fDz58+enYsWKdOrUyeuYItlGqs74mNkrwIPAzsCk55xz3wXm9QG6AgnAE865SYHpDYBYIDfwHdDDpSBEKJzx+eqrr5g4cSKxsbHkzJnT0ywi4nf8+HHi4uJ0N/U0FhkZSVRUFBEREcmm64zP2YXCcUpC39GjRxkzZgyffPIJM2bMwMy4/vrruffee7n99tvJnz+/1xFFMrX0PuPzjnOu3ykfWBPoANQCSgM/mllV51wCMBDoBszFX/jcCHyfBjnS3c6dO1mzZo2KHpEQEhcXR/78+alQoYLOxqYR5xy7du0iLi6OihUreh1HJEvJlSsXHTt2pGPHjqxfv54vv/ySoUOH0qVLFx555BFuv/12OnXqRNOmTQkPD/c6rkiWkl6DG7QGhjvnjjrnNgDrgCvMrBRQwDn3S+Asz1CgTTplSHNPPvkk+jZPJLQcOXKEokWLquhJQ2ZG0aJFdRZNJJ1deumlvPTSS6xdu5bZs2fTqVMnJkyYwA033EDp0qXp3r0706ZN06hwImkkLQqfx8zsNzMbYmaFA9PKAJuSLBMXmFYm8PzU6SIiF01FT9rTNhXJOGZG48aN+fjjj9m6dSujRo3C5/MxdOhQmjRpQunSpXnkkUdUBImk0nkLHzP70cyWn+HRGn+3tUpAXWArkDgsxpmOmO4c08/22d3MbIGZLdi5c+fZFssQO3bsoGHDhkyePNnTHCIiIpJ1RUZG0rZtW0aMGMHOnTsZOXIkPp+Pzz//nCZNmlCqVCm6du3Kt99+y8GDB72OK5KpnPcaH+dcs5S8kZkNBhLH44sDyiaZHQVsCUyPOsP0s332IGAQ+C8aTUmO9LJv3z7y5MlDZGSklzFEREQkm8iTJw/t2rWjXbt2HDx4kO+//57Ro0czevRohgwZQmRkJM2aNaNVq1a0bNmSUqVKeR1ZJKSlqqtb4JqdRLcBywPPxwEdzCyXmVUEqgDznHNbgf1m1tD8/SjuBb5NTYaMUrlyZX766SeuvfZar6OISAh64403qFWrFnXq1KFu3br8+uuvAISHh1O3bt3go2/fvgAcPnyY6667Llm3lZdeeono6Ohk90w4l2PHjnHttddy4sSJFOfctGkT119/PTVq1KBWrVq89957F7imIuKFvHnz0q5dO4YNG8bOnTuZOnUqDz30EMuXL6dbt26ULl2ayy+/nBdeeIFZs2Zd0H5BJLtI7ahub5tZXfzd1f4EHgJwzv1uZt8AK4ATwKOBEd0AuvO/4ay/J5OM6CYicja//PILEyZMYNGiReTKlYv4+HiOHTsGQO7cuVmyZMlprxkyZAi33357cNSmSZMmsXjxYpYsWcKKFSt49tln6dat2zk/N2fOnDRt2pQRI0Zw9913pyhrjhw56N+/P/Xr12f//v00aNCA5s2bU7NmzQtbaRHxTEREBE2aNKFJkya88847LF++nHHjxvH999/Tt29f3njjDQoWLEjTpk258cYbueGGGyhXrpzXsUU8l6ozPs65Ts65aOdcHedcq8AZncR5bzjnKjnnqjnnvk8yfYFzrnZg3mMpuYdPKLjrrrvO+0eIiGRPW7dupVixYuTKlQuAYsWKUbp06XO+5quvvqJ169bB9rhx4+jSpQvHjx/ngw8+oG3btin67DZt2vDVV1+dNn3+/PnUqVOHI0eOcPDgQWrVqsXy5cspVaoU9evXByB//vzUqFGDzZs3p3RVRSTEmBnR0dE8//zzzJo1i/j4eEaNGkX79u2ZN28e3bp1o3z58lSpUoWHHnqIESNGsH37dq9ji3givYazznIqVKigb0tEMgGfz0dsbCzgv7mpz+fjyy+/BODQoUP4fD5GjBgBwN69e/H5fIwZMwaA+Ph4fD4f48ePB2Dbtm0p+swWLVqwadMmqlatyiOPPMKMGTOC8w4fPpysq9uIESM4duwY69evp0KFCsHlFi5cyP79+ylatCizZs2iY8eOKfrs2rVrM3/+/NOmX3755bRq1YoXXniBZ555hnvuuYfatWsnW+bPP/9k8eLFXHnllSn6LBEJfYUKFaJt27YMHjyYjRs38vvvvzNgwABq1KjB8OHD6dChA5dccgm1a9fmiSeeYMyYMezYscPr2CIZQoVPCr355pu88MILXscQkRCUL18+Fi5cyKBBgyhevDh33nlnsPhK7OqW+LjzzjuJj4+nUKFCwdefPHmSuLg4unTpQnx8PA0aNGDAgAEA3HbbbQAMHTo0WJCdOHGCW265hSNHjhAeHk7OnDnZv3//ableeuklpkyZwoIFC3jmmWeSzTtw4ABt27bl3XffpUCBAumwVUTEa2ZGzZo1eeqppxg3bhy7du3i119/pW/fvpQpU4ZPP/2Utm3bUrJkSapVq0bXrl2JjY1l3bp1ZJIOOSIXJLXX+GQLzjnd00Ikk5g+fXrweURERLJ2njx5krULFiyYrF2sWLFk7UsuuSTFnxseHo7P58Pn8xEdHc3nn39Oly5dzrhs7ty5k90cdPXq1VSpUiU476qrrmLbtm1s2rSJMmX8tzqrU6cO3333HbfeeisDBw7kwQcfDI4yefTo0TOOOLl7924OHDjA8ePHOXLkCHnz5gX8Z8Latm3L3Xffze23357idRSRzC1HjhxcccUVXHHFFTz77LMcO3aMhQsX8vPPPzNr1izGjh3LkCFDAP/+r1GjRlx55ZU0bNiQmJiY4D5EJLPSGZ8U+Oqrr4iKimLjxo1eRxGRELR69WrWrl0bbC9ZsoTy5cufdfnChQuTkJAQLH4WL17M0aNHSUhI4OjRo3z99de0adOGhQsX0qBBAwBq1qzJqlWr2L17N3PmzKFNmzYA7Nq1i+LFixMREXHa53Tr1o3XXnuNu+++m2effRbwf5HTtWtXatSoQc+ePdNqE4hIJpQzZ04aNWrEM888w7hx44iPj2f58uV8/PHHNG3alGXLltG7d298Ph8FCxakXr16dO/endjYWJYvX66R4yTT0RmfFChbtixNmza9oG9/RST7OHDgAI8//jh79uwhR44cVK5cOTgcdeI1PoluvPFG+vbtS4sWLZg1axbNmjVjyZIlHD58mEqVKlGsWDEeeeQRLrvsMkaNGkW7du0A/x8oR48e5Y033uDll18Ovt+0adO4+eabT8s0dOhQcuTIwV133UVCQgKNGzfmp59+ImfOnHzxxRdER0cHc7355ptnfA8RyV7CwsKoVasWtWrV4qGHHgJg586dzJs3j7lz5zJ37ly+/vprPv74Y8B/hrpevXrExMTQoEEDYmJiqFatWnC0SpFQY5mlD2dMTIxbsGCB1zFEJMSsXLmSGjVqeB3jgi1evJgBAwbwxRdf0Lx5c955553TBh+4/fbb+eabb8iRw/8d1b333kuZMmV46623ki3z1ltvUa1atTTPeKZta2YLnXMxaf5hWYCOU5IdnDx5ktWrV7Nw4UIWLFjAwoULWbRoEYcOHQL8xVDiFyuJj+joaPLly+dxcskuznWc0hkfEREP1KtXj+uvv56EhARWrVpF9erVT1smcbS5REOHDk3WPnbsGG3atEmXokdE5EzCwsKoUaMGNWrU4J577gEI7scWLFgQHMhl5MiRwTPfZkblypWJjo6mdu3a1K5dm+joaCpXrhz8YkckI+iMTwo0btyYsmXLBofAFZHQkVnP+GQGOuNzYXTGR+R/nHNs2rSJJUuWsHTpUpYsWcLy5ctZt24dJ0+eBCBXrlzUqFGDWrVqBYupmjVrUqlSpTNetyiSEjrjk0rt27enYMGCXscQERERyRTMjHLlylGuXDlatWoVnH748GFWrlzJ8uXLWbZsGcuWLWPmzJnJbsQcERFB5cqVqVmzJtWqVaNq1arBn0WKFPFidSSLUOGTAk899ZTXEUREREQyvdy5c1O/fn3q16+fbPqBAwdYtWoVK1euZMWKFaxcuZJly5bx7bffJhs9rlixYlStWpWqVatSuXLl4KNSpUrJ7o8mciYqfM7DOcfJkyc1QomIiIhIOsmXLx8xMTHExCTvoXT8+HE2bNjAmjVrWL16dfDnpEmTgjeKTlSsWDEqVapEpUqVqFixIhUrVuTSSy+lYsWKREVF6XoiUeFzPnv37qVIkSK8//77PPbYY17HEREREck2IiIigmd4WrZsmWzewYMHWb9+PX/88Qfr1q0LPn755RdGjBhBQkJCcNkcOXJQrlw5ypcvT4UKFU77WaZMGV1XlA2o8DmPsLAwXnzxRS6//HKvo4iIiIhIQN68eYmOjiY6Ovq0eSdOnGDTpk1s2LCB9evXB3/+9ddfTJo0iS1btiRbPiwsjEsuuYRy5cpRtmxZypYtG3weFRVFVFQUJUuWVA+gTE6Fz3kUKFCAV1991esYIiIiIpJCOXLkCHZ3a9KkyWnzjx49yqZNm/jrr7/4888/2bhxI5s2bWLjxo0sXbqU8ePHc+TIkWSvCQ8Pp1SpUpQpU4aoqCjKlClD6dKlkz1KlSpFwYIFMbOMWlW5ACp8zuP48eMAOv0pImeVL18+Dhw4cMZ5PXr0YNSoUWzatInff/+dTp06AbBx40YKFixIwYIFKVasGJ9++ik1atRIdk+enj17cu+992bIOoiIZCe5cuUKDoxwJs45du3axcaNG9m8eTObN28mLi4u+HPFihVMmTKFffv2nfba3LlzU6pUKS655JLTfl5yySWULFmSkiVLUqJECXLlypXeqypJqPA5j7Fjx3LnnXeyfPlyatWq5XUcEclETp48ydixYylbtiwzZ87E5/OxZMkSALp06ULLli1p164dAH/++SeVKlUKzhcREe+YGcWKFaNYsWKnjUCX1MGDB9m6dStbtmwJPjZv3sy2bdvYtm0bK1as4KeffuLvv/8+4+sLFSqUrBBK+ihevHjwZ/HixSlcuLC62qWSCp/zqFWrFv/85z8pXbq011FEJJOZNm0atWvX5s4772TYsGH4fD6vI4mISBrKmzfvOc8cJTpy5Ajbt29n27ZtbN++/YyP33//nWnTprFr164zvkdYWBhFihQJFmTFixenWLFiFC1aNPjz1OeFChVSsZSECp/zqFWrls70iGQiKSkuWrZsSa9evYLLd+nShS5duhAfHx88A5No+vTpF51l2LBhdOzYkdatW/Pcc89x/Pjxc3ab/eOPP6hbt26w/Z///Idrrrnmoj9fRERCQ2RkJOXLl6d8+fLnXfbEiRPEx8ezc+dOduzYwY4dO4LtpD/XrFnDnDlz2LVrV7J7HZ2qUKFCFC1alCJFigQfhQsXDv5MfCS2CxUqROHChcmbN2+Wu1ZJhc957N+/n4iICCIjI72OIiKZyLFjx/juu+945513yJ8/P1deeSWTJ0/mlltuOetr1NVNRERy5MgRvB4oJZxz7N+/n127dhEfH8+uXbvYtWsXu3fvDv5M+nzdunX8/fff7Nmzh5MnT571fcPDw4NFUKFChShUqBAFCxY87Xni9aoFChQIPk9s58yZM602S5pQ4XMePXv2ZOLEiacNeygioelCz9AkXb5YsWKpOsOT1A8//MDevXuDw6weOnSIPHnynLPwERERuVBmRoECBShQoAAVK1ZM8etOnjzJ/v372b17N3///XfwsWfPHvbs2RN8njh97969bNmyhT179rB3714OHjx43s+IjIwMZjv1kT9//uDPpM9r166dbKCftKTC5zzat2+ve/iIyAUbNmwYn376KR07dgT8F8BWrFgxWACJiIh4KSwsLHh25kIKpkTHjx9n7969Z3zs27cv+DPpY+/evfz555/s27eP/fv3s2/fvuAIyolefPFF/vnPf6bVaiajwuc8WrRo4XUEEQlxhw4dIioqKth+5JFHmDRpEp988klwWt68ebn66qsZP348d9555xnf59RrfO6//36eeOKJdMstIiJysSIiIoIDLaTG0aNHg0XQ/v37U/1+56LC5zy2bdtGgQIF9A2tiJzVmfpIP/fcc6dNGzNmTPB5bGxssnkVKlTg8OHDaZ5NREQklOXKlYtcuXKla8GTKCzdPyGTu/LKK+nevbvXMUREREREJBV0xuc8Xn/9dcqUKeN1DBERERERSQUVPufRqVMnryOIiIiIiEgqqavbORw/fpz169er371IiHPOeR0hy9E2FRGRrEaFzzmsX7+eSpUqMXbsWK+jiMhZREZGsmvXLv2hnoacc+zatSvb37jZzG40s9Vmts7MenudR0REUkdd3c6hRIkSxMbG0qhRI6+jiMhZREVFERcXx86dO72OkqVERkYmG6I7uzGzcOBDoDkQB8w3s3HOuRXeJhMRkYulwuccChcuTOfOnb2OISLnEBERcVE3XhM5jyuAdc659QBmNhxoDajwERHJpNTV7Rzi4+NZvXo1J06c8DqKiIhkrDLApiTtuMA0ERHJpFT4nMOwYcOoXr06e/bs8TqKiIhkLDvDtNMuJDOzbma2wMwWqLuliEhoU+FzDjfeeCNffvklhQoV8jqKiIhkrDigbJJ2FLDl1IWcc4OcczHOuZjixYtnWDgREblwlllGQjKzncBfqXiLYkB8GsXJjLL7+oO2AWgbZPf1h9Rvg/LOuSz/F76Z5QDWAE2BzcB84C7n3O/neE1WP06Fej4I/YzKl3qhnlH5Ui/djlOZZnCD1B5ozWyBcy4mrfJkNtl9/UHbALQNsvv6g7ZBSjnnTpjZY8AkIBwYcq6iJ/CaLH2cCvV8EPoZlS/1Qj2j8qVeembMNIWPiIhIRnLOfQd853UOERFJG7rGR0REREREsrzsVPgM8jqAx7L7+oO2AWgbZPf1B22DUBbq/zahng9CP6PypV6oZ1S+1Eu3jJlmcAMREREREZGLlZ3O+IiIiIiISDaVpQofM7vRzFab2Toz632G+WZm7wfm/2Zm9b3ImZ5SsA3uDqz7b2Y2x8wu8yJnejrfNkiy3OVmlmBm7TIyX3pLyfqbmc/MlpjZ72Y2I6MzprcU/D8oaGbjzWxpYBvc50XO9GJmQ8xsh5ktP8v8LL8vDGWhvp9OQb7WgWxLAjdvvTqU8iVZzrN9fAq2oc/M9ga24RIzeymU8iXJ6MlxIgXb7x9Jtt3ywL9zkRDL6OlxJgX5CpvZ2MD/5XlmVjuD83lznHLOZYkH/uFG/wAuBXICS4GapyxzM/A9/jtyNwR+9Tq3B9ugMVA48Pym7LgNkiz3E/4Rm9p5nTuDfwcKASuAcoF2Ca9ze7ANngP+FXheHNgN5PQ6expug2uB+sDys8zP0vvCUH6E+n46hfny8b+u8nWAVaGUL8lynuzjU7gNfcCEEP4d9Ow4kdJ/4yTL3wr8FILb0LPjTArz/Rt4OfC8OjA1g7ehJ8eprHTG5wpgnXNuvXPuGDAcaH3KMq2Boc5vLlDIzEpldNB0dN5t4Jyb45z7O9Cci/9u5FlJSn4PAB4HRgM7MjJcBkjJ+t8FjHHObQRwzmXHbeCA/GZm+P+I2w2cyNiY6cc5NxP/Op1NVt8XhrJQ30+nJN8BF/jLBMiL//9TyOQL8HIfn9KMXgn148SFbr+OwLAMSfY/oX6cSUm+msBUAOfcKqCCmZXMoHyeHaeyUuFTBtiUpB0XmHahy2RmF7p+XfFX01nJebeBmZUBbgM+zsBcGSUlvwNVgcJmNt3MFprZvRmWLmOkZBt8ANQAtgDLgB7OuZMZEy8kZPV9YSgL9f10ivKZ2W1mtgqYCNyfQdkgc+zjU/pv3CjQDep7M6uVMdGA0D9OpPj/iJnlAW7EX+RmpFA/zqQk31LgdgAzuwIoT2h9GZ4ux6msdANTO8O0U7+FSskymVmK18/Mrsd/QM3QvtkZICXb4F3gWedcgv+LmCwlJeufA2gANAVyA7+Y2Vzn3Jr0DpdBUrINbgCWAE2ASsAUM/vZObcvnbOFiqy+Lwxlob6fTlE+59xYYKyZXQu8BjRL72ABmWEfn5KMi4DyzrkDZnYz8H9AlfQOFhDqx4kL2T/dCsx2zp3rzEF6CPXjTEry9QXeM7Ml+AuzxYRWz4d0OU5lpcInDiibpB2Fv8q+0GUysxStn5nVAT4FbnLO7cqgbBklJdsgBhgeOCAWA242sxPOuf/LkITpK6X/D+KdcweBg2Y2E7gMyCqFT0q2wX1A30B3nXVmtgF/H+d5GRPRc1l9XxjKQn0/fUG/G865mWZWycyKOefi0z1d5tjHnzdj0j9+nXPfmdlHIbYNvTxOXMjvYAcyvpsbhP5xJqW/g/eBfyABYEPgESrS5ziVHhcsefHAX8StByryvwu5ap2yzC0kv1Bqnte5PdgG5YB1QGOv83q1DU5ZPpasNbhBSn4HauDv15sDyAMsB2p7nT2Dt8FA4JXA85LAZqCY19nTeDtU4OwXjWbpfWEoP0J9P53CfJX53+AG9QP/fyxU8p2yfIbv41O4DS9Jsg2vADaG0jb08jiR0n9joCD+a0TyZuS/7wVsQ8+OMynMV4jAYAvAg/ivp8no7Zjhx6ksc8bHOXfCzB4DJuEfzWKIc+53M3s4MP9j/KO73Iz/gHKIQKWbVaRwG7wEFAU+CnwbdsI5F+NV5rSWwm2QZaVk/Z1zK83sB+A34CTwqXPujMNJZkYp/B14DYg1s2X4d6rPuoz5pjVDmNkw/KNGFTOzOOBlIAKyx74wlIX6fjqF+doC95rZceAwcKcL/KUSIvk8lcKM7YDuZnYC/zbsEErb0MvjxAX8G98GTHb+s1IZKtSPMynMVwMYamYJ+Efw65oR2RJ5dZyyDPp/JiIiIiIi4pmsNKqbiIiIiIjIGanwERERERGRLE+Fj4iIiIiIZHkqfEREREREJMtT4SMiIiIiIlmeCh/JtszsEjMbbmZ/mNkKM/vOzKqm02f5zGxC4HkrM+udHp8jIiJZi45VImkny9zHR+RCBO5SPBb43DnXITCtLv6bjKXrnamdc+OAcen5GSIikvnpWCWStnTGR7Kr64HjSW9255xbAswys3+b2XIzW2Zmd0LwW7AZZvaNma0xs75mdreZzQssVymwXKyZfWxmPweWa3nqB5tZFzP7IPD8VjP71cwWm9mPZlYyMP0VMxtiZtPNbL2ZPZHk9fea2W9mttTMvghMK25mo81sfuBxVTpuOxERyRg6VomkIZ3xkeyqNrDwDNNvB+oClwHFgPlmNjMw7zL8dzreDazHfyfrK8ysB/A48GRguQrAdUAlYJqZVT5HjllAQ+ecM7MHgGeApwPzquM/6OUHVpvZQKAq8DxwlXMu3syKBJZ9D3jHOTfLzMrhv1tzjRRuCxERCU06VomkIRU+IsldDQxzziUA281sBnA5sA+Y75zbCmBmfwCTA69Zhn+nn+gb59xJYK2Zrcd/UDibKGCEmZUCcgIbksyb6Jw7Chw1sx34uzY0AUY55+IBnHO7A8s2A2r6e0UAUMDM8jvn9l/4JhARkRCnY5XIRVBXN8mufgcanGG6nWFaoqNJnp9M0j5J8i8R3CmvO7Wd1H+AD5xz0cBDQORZPi8h8Bl2lvcLAxo55+oGHmV0IBERyfR0rBJJQyp8JLv6CchlZg8mTjCzy4G/gTvNLNzMigPXAvMu8L3bm1lYoC/1pcDqcyxbENgceN45Be89FbjDzIoGMid2H5gMPJa4UODiVxERydx0rBJJQyp8JFtyzjngNqC5+YcI/R14Bfga+A1Yiv+A84xzbtsFvv1qYAbwPfCwc+7IOZZ9BRhpZj8D8SnI/TvwBjDDzJYCAwKzngBiAheSrgAevsDMIiISYnSsEklb5v8/JSJpwcxigQnOuVFeZxERETkTHasku9IZHxERERERyfJ0xkdERERERLI8nfEREREREZEsT4WPiIiIiIhkeSp8REREREQky1PhIyIiIiIiWZ4KHxERERERyfJU+IiIiIiISJb3/5LHNg5ybdu5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1008x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "se_formula_iv = lambda compliance: np.std(e_iv)/(compliance*np.std(Z)*np.sqrt(n))\n",
    "x = np.linspace(0.01, 1, 50)\n",
    "\n",
    "effect = iv_regr.params[\"prime_card\"]\n",
    "\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))\n",
    "\n",
    "ax1.plot(x, effect-se_formula_iv(x)*2, label=\"SE($\\\\beta_{IV}$) x2\", ls=\":\", color=\"0\")\n",
    "ax1.plot(x, effect+se_formula_iv(x)*2, ls=\":\", color=\"0\")\n",
    "ax1.hlines(effect, 0, 1, ls=\"-.\", label=\"LATE\")\n",
    "ax1.hlines(0, 0, 1)\n",
    "ax1.set_xlabel(\"Compliance\")\n",
    "ax1.set_ylim(-(effect+100), (effect+100)*2)\n",
    "ax1.legend(loc=\"lower right\")\n",
    "ax1.set_title(\"95% CI around LATE\");\n",
    "\n",
    "\n",
    "x = np.linspace(0.2, 1, 50)\n",
    "ax2.plot(x, 1/(x**2))\n",
    "ax2.hlines(10, 0.2, 1, ls=\":\", label=\"10\")\n",
    "ax2.hlines(4, 0.2, 1, ls=\"-.\", label=\"4\")\n",
    "ax2.set_xlabel(\"Compliance\")\n",
    "ax2.set_ylabel(\"$N_{IV}$/N\")\n",
    "ax2.set_title(\"Required Sample Size Ratio\")\n",
    "ax2.legend()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b15a0ab",
   "metadata": {},
   "source": [
    "## Additional Controls and Instruments\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "a520339f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:14.191876Z",
     "start_time": "2024-01-23T12:22:14.018689Z"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.50.0 (0)\n",
       " -->\n",
       "<!-- Pages: 1 -->\n",
       "<svg width=\"270pt\" height=\"206pt\"\n",
       " viewBox=\"0.00 0.00 270.49 206.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 202)\">\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-202 266.49,-202 266.49,4 -4,4\"/>\n",
       "<!-- U -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>U</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"30.55\" cy=\"-180\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"30.55\" y=\"-176.3\" font-family=\"Times,serif\" font-size=\"14.00\">U</text>\n",
       "</g>\n",
       "<!-- T -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>T</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"134.79\" cy=\"-124\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"134.79\" y=\"-120.3\" font-family=\"Times,serif\" font-size=\"14.00\">T</text>\n",
       "</g>\n",
       "<!-- U&#45;&gt;T -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>U&#45;&gt;T</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M52.24,-168.67C67.3,-160.42 87.91,-149.13 104.64,-139.96\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"106.39,-143 113.48,-135.12 103.03,-136.86 106.39,-143\"/>\n",
       "</g>\n",
       "<!-- Y -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>Y</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"235.49\" cy=\"-97\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"235.49\" y=\"-93.3\" font-family=\"Times,serif\" font-size=\"14.00\">Y</text>\n",
       "</g>\n",
       "<!-- U&#45;&gt;Y -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>U&#45;&gt;Y</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M57.51,-178.61C86.81,-176.17 135.21,-169.46 172.49,-151 187.87,-143.38 202.48,-131.01 213.65,-120.02\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"216.31,-122.31 220.8,-112.72 211.31,-117.42 216.31,-122.31\"/>\n",
       "</g>\n",
       "<!-- T&#45;&gt;Y -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>T&#45;&gt;Y</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M159.94,-117.4C172.12,-114.06 187.11,-109.96 200.45,-106.31\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"201.53,-109.65 210.25,-103.63 199.68,-102.89 201.53,-109.65\"/>\n",
       "</g>\n",
       "<!-- Z -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>Z</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"30.55\" cy=\"-126\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"30.55\" y=\"-122.3\" font-family=\"Times,serif\" font-size=\"14.00\">Z</text>\n",
       "</g>\n",
       "<!-- Z&#45;&gt;T -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>Z&#45;&gt;T</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M57.61,-125.49C69.75,-125.25 84.39,-124.97 97.58,-124.71\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"97.84,-128.2 107.77,-124.51 97.71,-121.21 97.84,-128.2\"/>\n",
       "</g>\n",
       "<!-- Income -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>Income</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"134.79\" cy=\"-70\" rx=\"37.89\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"134.79\" y=\"-66.3\" font-family=\"Times,serif\" font-size=\"14.00\">Income</text>\n",
       "</g>\n",
       "<!-- Income&#45;&gt;Y -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>Income&#45;&gt;Y</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M167.84,-78.77C178.24,-81.61 189.83,-84.78 200.38,-87.67\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"199.54,-91.07 210.11,-90.33 201.39,-84.32 199.54,-91.07\"/>\n",
       "</g>\n",
       "<!-- Age -->\n",
       "<g id=\"node6\" class=\"node\">\n",
       "<title>Age</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"30.55\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"30.55\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">Age</text>\n",
       "</g>\n",
       "<!-- Age&#45;&gt;T -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>Age&#45;&gt;T</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M48.13,-31.95C52.54,-35.98 57.16,-40.5 61.09,-45 79.59,-66.17 77.51,-76.84 97.09,-97 99.93,-99.92 103.1,-102.75 106.37,-105.42\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"104.46,-108.36 114.53,-111.65 108.71,-102.8 104.46,-108.36\"/>\n",
       "</g>\n",
       "<!-- Age&#45;&gt;Y -->\n",
       "<g id=\"edge7\" class=\"edge\">\n",
       "<title>Age&#45;&gt;Y</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M57.76,-18.5C87.06,-19.98 135.27,-25.2 172.49,-43 188.18,-50.5 202.99,-63.1 214.19,-74.24\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"211.87,-76.88 221.34,-81.64 216.91,-72.02 211.87,-76.88\"/>\n",
       "</g>\n",
       "<!-- Score -->\n",
       "<g id=\"node7\" class=\"node\">\n",
       "<title>Score</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"30.55\" cy=\"-72\" rx=\"30.59\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"30.55\" y=\"-68.3\" font-family=\"Times,serif\" font-size=\"14.00\">Score</text>\n",
       "</g>\n",
       "<!-- Score&#45;&gt;T -->\n",
       "<g id=\"edge8\" class=\"edge\">\n",
       "<title>Score&#45;&gt;T</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M54.24,-83.54C68.76,-90.92 87.76,-100.59 103.55,-108.62\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"102.4,-111.96 112.9,-113.38 105.58,-105.72 102.4,-111.96\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.graphs.Digraph at 0x7f84aa427810>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gr = Digraph(format=\"png\", graph_attr={\"rankdir\":\"LR\"})\n",
    "\n",
    "gr.edge(\"U\", \"T\")\n",
    "gr.edge(\"U\", \"Y\")\n",
    "gr.edge(\"Z\", \"T\")\n",
    "gr.edge(\"T\", \"Y\")\n",
    "gr.edge(\"Income\", \"Y\")\n",
    "gr.edge(\"Age\", \"T\")\n",
    "gr.edge(\"Age\", \"Y\")\n",
    "gr.edge(\"Score\", \"T\")\n",
    "\n",
    "gr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "304fde6b",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:14.291008Z",
     "start_time": "2024-01-23T12:22:14.193566Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Parameter Estimates</caption>\n",
       "<tr>\n",
       "       <td></td>      <th>Parameter</th> <th>Std. Err.</th> <th>T-stat</th> <th>P-value</th> <th>Lower CI</th> <th>Upper CI</th>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>   <td>2519.4</td>    <td>21.168</td>   <td>119.02</td> <td>0.0000</td>   <td>2477.9</td>   <td>2560.9</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>prime_card</th>  <td>659.04</td>    <td>58.089</td>   <td>11.345</td> <td>0.0000</td>   <td>545.19</td>   <td>772.90</td> \n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula = 'pv ~ 1 + [prime_card ~ prime_elegible + credit_score]'\n",
    "iv_model = IV2SLS.from_formula(formula, df).fit()\n",
    "\n",
    "iv_model.summary.tables[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "798c4e04",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:14.343456Z",
     "start_time": "2024-01-23T12:22:14.292449Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Parameter Estimates</caption>\n",
       "<tr>\n",
       "       <td></td>      <th>Parameter</th> <th>Std. Err.</th> <th>T-stat</th> <th>P-value</th> <th>Lower CI</th> <th>Upper CI</th>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>   <td>210.62</td>    <td>37.605</td>   <td>5.6008</td> <td>0.0000</td>   <td>136.91</td>   <td>284.32</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>age</th>         <td>9.7444</td>    <td>0.8873</td>   <td>10.982</td> <td>0.0000</td>   <td>8.0053</td>   <td>11.483</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>income</th>      <td>0.3998</td>    <td>0.0008</td>   <td>471.04</td> <td>0.0000</td>   <td>0.3981</td>   <td>0.4014</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>prime_card</th>  <td>693.12</td>    <td>12.165</td>   <td>56.978</td> <td>0.0000</td>   <td>669.28</td>   <td>716.96</td> \n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula = '''pv ~ 1 \n",
    "+ [prime_card ~ prime_elegible + credit_score]\n",
    "+ income + age'''\n",
    "\n",
    "iv_model = IV2SLS.from_formula(formula, df).fit(cov_type=\"unadjusted\")\n",
    "\n",
    "iv_model.summary.tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7da0f4c8",
   "metadata": {},
   "source": [
    "### 2SLS by Hand\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "4a3ca2f9",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:14.376607Z",
     "start_time": "2024-01-23T12:22:14.344543Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "       <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>  <td>  210.6177</td> <td>   40.832</td> <td>    5.158</td> <td> 0.000</td> <td>  130.578</td> <td>  290.657</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>prime_card</th> <td>  693.1207</td> <td>   13.209</td> <td>   52.474</td> <td> 0.000</td> <td>  667.229</td> <td>  719.013</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>income</th>     <td>    0.3998</td> <td>    0.001</td> <td>  433.806</td> <td> 0.000</td> <td>    0.398</td> <td>    0.402</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>age</th>        <td>    9.7444</td> <td>    0.963</td> <td>   10.114</td> <td> 0.000</td> <td>    7.856</td> <td>   11.633</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula_1st = \"prime_card ~ prime_elegible + credit_score + income+age\"\n",
    "first_stage = smf.ols(formula_1st, data=df).fit()\n",
    "\n",
    "iv_model = smf.ols(\n",
    "    \"pv ~ prime_card + income + age\",\n",
    "    data=df.assign(prime_card=first_stage.fittedvalues)).fit()\n",
    "\n",
    "iv_model.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73042769",
   "metadata": {},
   "source": [
    "### Matrix Implementation\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "3d789d29",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:14.389797Z",
     "start_time": "2024-01-23T12:22:14.378246Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([693.12072518])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Z = df[[\"prime_elegible\", \"credit_score\", \"income\", \"age\"]].values\n",
    "X = df[[\"prime_card\", \"income\", \"age\"]].values\n",
    "Y = df[[\"pv\"]].values\n",
    "\n",
    "def add_intercept(x):\n",
    "    return np.concatenate([np.ones((x.shape[0], 1)), x], axis=1)\n",
    "\n",
    "Z_ = add_intercept(Z)\n",
    "X_ = add_intercept(X)\n",
    "\n",
    "# pre-multiplying Z_.dot(...) last is important to avoid\n",
    "# creating a huge NxN matrix\n",
    "X_hat = Z_.dot(np.linalg.inv(Z_.T.dot(Z_)).dot(Z_.T).dot(X_))\n",
    "\n",
    "b_iv = np.linalg.inv(X_hat.T.dot(X_hat)).dot(X_hat.T).dot(Y)\n",
    "b_iv[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "d50d8d46",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:14.395998Z",
     "start_time": "2024-01-23T12:22:14.391505Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "12.164694395033125"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e_hat_iv = (Y - X_.dot(b_iv))\n",
    "\n",
    "var = e_hat_iv.var()*np.diag(np.linalg.inv(X_hat.T.dot(X_hat)))\n",
    "\n",
    "np.sqrt(var[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "9f4f9fc3",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:14.410300Z",
     "start_time": "2024-01-23T12:22:14.397376Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "12.156252763192523"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_tilde = smf.ols(\"prime_card ~ income + age\", data=df).fit().resid\n",
    "\n",
    "e_hat_iv.std()/(t_tilde.std()*np.sqrt(n*first_stage.rsquared))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0f3a4cc",
   "metadata": {},
   "source": [
    "## Discontinuity Design\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d99e943",
   "metadata": {},
   "source": [
    "### Discontinuity Design Assumptions\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2177b348",
   "metadata": {},
   "source": [
    "### Intention to Treat Effect\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "2fd1046d",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:14.428637Z",
     "start_time": "2024-01-23T12:22:14.411988Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>balance</th>\n",
       "      <th>prime_card</th>\n",
       "      <th>pv</th>\n",
       "      <th>tau</th>\n",
       "      <th>categ</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>12100.0</td>\n",
       "      <td>1</td>\n",
       "      <td>356.472</td>\n",
       "      <td>300.0</td>\n",
       "      <td>always-takers</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4400.0</td>\n",
       "      <td>1</td>\n",
       "      <td>268.172</td>\n",
       "      <td>300.0</td>\n",
       "      <td>always-takers</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4600.0</td>\n",
       "      <td>1</td>\n",
       "      <td>668.896</td>\n",
       "      <td>300.0</td>\n",
       "      <td>always-takers</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3500.0</td>\n",
       "      <td>1</td>\n",
       "      <td>428.094</td>\n",
       "      <td>300.0</td>\n",
       "      <td>always-takers</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>12700.0</td>\n",
       "      <td>1</td>\n",
       "      <td>1619.793</td>\n",
       "      <td>700.0</td>\n",
       "      <td>complier</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   balance  prime_card        pv    tau          categ\n",
       "0  12100.0           1   356.472  300.0  always-takers\n",
       "1   4400.0           1   268.172  300.0  always-takers\n",
       "2   4600.0           1   668.896  300.0  always-takers\n",
       "3   3500.0           1   428.094  300.0  always-takers\n",
       "4  12700.0           1  1619.793  700.0       complier"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_dd = pd.read_csv(\"./data/prime_card_discontinuity.csv\")\n",
    "df_dd.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "a1c54382",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:14.453943Z",
     "start_time": "2024-01-23T12:22:14.430325Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "                 <td></td>                   <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>                      <td>  559.2980</td> <td>    8.395</td> <td>   66.621</td> <td> 0.000</td> <td>  542.843</td> <td>  575.753</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>I(balance > 0)[T.True]</th>         <td>  261.0699</td> <td>   10.128</td> <td>   25.777</td> <td> 0.000</td> <td>  241.218</td> <td>  280.922</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>balance</th>                        <td>    0.0616</td> <td>    0.005</td> <td>   11.892</td> <td> 0.000</td> <td>    0.051</td> <td>    0.072</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>balance:I(balance > 0)[T.True]</th> <td>   -0.0187</td> <td>    0.005</td> <td>   -3.488</td> <td> 0.000</td> <td>   -0.029</td> <td>   -0.008</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = smf.ols(f\"pv~balance*I(balance>0)\",\n",
    "            df_dd.assign(balance = lambda d: d[\"balance\"] - 5000)).fit()\n",
    "m.summary().tables[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "eefb7795",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:14.636458Z",
     "start_time": "2024-01-23T12:22:14.455842Z"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f84bb70a890>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfsAAAEKCAYAAAAVRfxuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABqFklEQVR4nO3deXzU1b34/9eZNcts2ROygiYsKi4gCiibyObS1tbr1orWXq21t2qvXvXb+2i9/d0+emvtct3a2kXFtbZq9QqIoggCIoJSUaiyQwghZJ/JJLOe3x8z8+kMhBAge97Px4MHyeczy8mHIe/POed93kdprRFCCCHE0GXq7wYIIYQQondJsBdCCCGGOAn2QgghxBAnwV4IIYQY4iTYCyGEEEOcBHshhBBiiOu1YK+UKlVKrVBKbVVKfaaUuj1+/OdKqX8opT5RSr2ilPLEj1copdqVUpvif36b9FoTlFKblVLblVIPKaVUb7VbCCGEGGpUb62zV0oVAUVa64+UUk5gI/BloAR4R2sdVkr9DEBrfY9SqgJ4XWt9eievtR64HVgHLAEe0lov7ZWGCyGEEENMr/XstdYHtNYfxb/2AluBYq31m1rrcPxh64gF/6OK3zS4tNbv69idySJiNw1CCCGE6AZLX7xJvNd+NvDBYae+Cfw56fuRSqmPgVbgP7XW7wHFQHXSY6rjx7qUm5urKyoqTqLVQgghxOCxcePGeq11Xmfnej3YK6UcwEvAHVrr1qTjPwDCwLPxQweAMq11g1JqAvA3pdRpQGfz853OPSilbgZuBigrK2PDhg0994MIIYQQA5hSas/RzvVqNr5Sykos0D+rtX456fhC4FLguvjQPFrrgNa6If71RmAHUEWsJ5881F8C1HT2flrrx7XWE7XWE/PyOr25EUIIIYad3szGV8Afga1a618mHZ8H3ANcrrX2Jx3PU0qZ41+PAiqBnVrrA4BXKXV+/DWvB17trXYLIYQQQ01vDuNPBb4BbFZKbYof+3/AQ4AdeCu+gm6d1vrbwDTgx0qpMBABvq21bow/71bgSSAdWBr/I4QQQohu6LVgr7VeTefz7UuO8viXiA35d3ZuA3DEkjwhhBBCHJtU0BNCCCGGuD5ZejfQRKNRqquraWtr6++miF6UmZlJSUkJJpPc0wohhrdhGezr6+tRSjF69GgJBENUNBpl//791NfXk5+f39/NEWLYi0ajNDY24vf7cTgcZGVlIZXP+86wjHTNzc0UFBRIoB/CTCYTBQUFtLS09HdThBBAS0sLe/bs4dChQ+zatQuv19vfTRpWhmW0i0QiWK3W/m6G6GVWq5VwOHzsBwohet3h/xdDoVA/tWR4GpbBHpDho2FA/o2FGDgcDgd2ux2AjIwMHA5HP7doeBmWc/ZCCCH6Vnp6OlVVVQSDQex2u4yu9jEJ9kIIIfqEzWbDZrP1dzOGpWE7jC+EEKJvaK3p6OggEAj0d1OGLQn2g8iBAwdYuHAheXl5pKWlMW7cOFauXGmc11pz//33M2LECNLT05kxYwafffZZP7ZYCDHcaa2pr6/ns88+47PPPpMVMv1Egv0g0dzczNSpU9Fas3jxYrZu3crDDz+csob8gQce4Be/+AUPP/wwH374Ifn5+Vx88cWyxEUI0W9CoRD79u0DYoG/pqaG+Ganog/JnP0J8nq9vPzyy2zevJkzzjiDK664AqfT2Wvv98ADD1BUVMSiRYuMYyNHjjS+1lrz61//mnvvvZevfvWrADz11FPk5+fz3HPPccstt/Ra24QQ4mhMJhMWi8VYamez2WSlTD+Qnv0J8Hq9XHPNNTz++OO8++67PP7441xzzTW92oP+29/+xnnnncdVV11Ffn4+Z511Fo888ohxh7xr1y5qa2uZM2eO8Zz09HSmTZvG2rVre61dQgjRFYvFwqhRo8jKyiI3N5cRI0b0d5OGJQn2J+Dll1+mqanJSDYJBAI0NTXx8ssv99p77ty5k8cee4xRo0axbNkybr/9du69914effRRAGprawEoKChIeV5BQYFxTggh+oPD4WDUqFGUl5eTnp7e380ZlmQY/wRs3rz5iKzSQCDAp59+2mvvGY1GmThxIj/96U8BOPvss9m2bRuPPvoo3/3ud43HHT48prWWITMhhBjmpGd/As444wyjElSC3W7n9NNP77X3LCoqYty4cSnHxo4dy969ewEoLCwEOKIXX1dXd0RvXwghxPAiwf4EXHHFFWRlZRkB3263k5WVxRVXXNFr7zl16lQ+//zzlGNffPEF5eXlQCxZr7CwkLfeess439HRwXvvvceUKVN6rV1CCCEGvl4L9kqpUqXUCqXUVqXUZ0qp2+PHs5VSbymltsX/zkp6zn1Kqe1Kqc+VUnOTjk9QSm2On3tI9fO4tNPp5Pnnn+fmm29m5syZ3HzzzTz//PO9mo1/5513sm7dOn7yk5+wfft2/vKXv/DQQw9x2223AbHh+zvuuIP/+Z//4eWXX+bTTz/lhhtuwOFwcO211/Zau4QQQgwCWute+QMUAefEv3YCXwDjgAeAe+PH7wV+Fv96HPB3wA6MBHYA5vi59cBkQAFLgfnHev8JEyboo9myZctRzw1kr7/+uh4/fry22+26srJS/+///q+ORqPG+Wg0qn/0ox/pwsJCbbfb9bRp0/TmzZv7scX9b7D+WwshxPECNuijxESl+6i4gVLqVeCR+J8ZWusDSqki4F2t9Wil1H0AWuufxh+/DLgf2A2s0FqPiR+/Jv78LheOT5w4UW/YsKHTc1u3bmXs2LE98nOJgU3+rYUQw4VSaqPWemJn5/pkzl4pVQGcDXwAFGitDwDE/06UgCsG9iU9rTp+rDj+9eHHhRBCCNENvR7slVIO4CXgDq11a1cP7eSY7uJ4Z+91s1Jqg1Jqw6FDh46/sUIIIcQQ1KvBXillJRbon9VaJyrOHIwP3xP/uy5+vBooTXp6CVATP17SyfEjaK0f11pP1FpPzMvL67kfRAghhBjEejMbXwF/BLZqrX+ZdOo1YGH864XAq0nHr1ZK2ZVSI4FKYH18qN+rlDo//prXJz1HCCGEEMfQmxX0pgLfADYrpTbFj/0/4H+AF5VSNwF7gSsBtNafKaVeBLYAYeA2rXUk/rxbgSeBdGLZ+Et7sd1CCDHsRaNRlFJSgXOI6LVgr7VeTefz7QAXHeU5PwF+0snxDUDvlacTQghhaGpqorq6GpvNRklJCZmZmf3dJHGSpIKeEEIIQyAQYNeuXQSDQXw+H3V1dcd+khjwJNgLIYRIkVx/JRqN9mNLRE+RYC+EEMJgt9sZOXIkZrOZtLQ02UhriJAtboUQQqTIzs7G5XKhlMJsNvd3c0QPkJ79ILJq1Souv/xyiouLUUrx5JNPGudCoRD33HMP48ePJzMzk6KiIq699lpjC9yEQCDAv/3bv5Gbm0tmZiaXX3451dXVCCFEMovFIoF+CJFgfxLC4TDNzc2Ew+E+eT+fz8fpp5/O//7v/5Kenp5yzu/389FHH/GDH/yAjz76iFdffZV9+/Yxb968lPbdcccdvPTSSzz//PO89957tLa2cumllxKJRA5/OyGEEEOEDOOfAK01f/rTn1i0aBGhUAir1cr111/PN7/5zV5dk7pgwQIWLFgAwA033JByzu12p+xlD/C73/2O0047ja1bt3LGGWfQ0tLCH//4R5544gkuvvhiAJ5++mnKy8tZvnw5c+fORQghxNAjPfsT8Kc//YknnniCtrY2gsEgbW1tPPHEEzzxxBP93bQUra2xrQiysrIA2LhxI6FQiDlz5hiPKS0tZezYsaxdu7Zf2iiEEKL3SbA/TuFwmEWLFtHR0ZFyvKOjg6eeeqrPhvSPJRgM8u///u9cdtlllJTEthaora3FbDaTm5ub8tiCggJqa2v7o5lCCCH6gAzjHyefz0coFOr0XCgUwufz4fF4+rZRhwmHw3z961+nubmZ11577ZiP11pLSUwhhBjCpGd/nBwOB1artdNzVqsVh8PRxy1KFQ6Hueaaa/jkk094++23ycnJMc4VFhYSiUSor69PeU5dXZ2spRVCiCFMgv1xslgsXH/99aSlpaUcT0tLY+HChVgs/TdYEgqFuOqqq/jkk09YsWIFhYWFKecnTJiA1WpNSeSrrq5m69atTJkypa+bK4QQoo/IMP4JSGTdP/XUU0Y2/sKFC7nxxht79X19Ph/bt28HYiUs9+7dy6ZNm8jOzmbEiBFceeWVfPjhh/zf//0fSiljHt7tdpOeno7b7eamm27i7rvvJj8/n5ycHL7//e8zfvx4Zs+e3attF0II0X9Ucg3koWTixIl6w4YNnZ7bunUrY8eOPen3CIfD+Hw+HA5Hn/To3333XWbOnHnE8YULF3L//fczcuTITp/3xBNPGEv1Ojo6uPvuu3nuuedob2/noosu4rHHHqO0tLQ3m95veurfWgghBjql1Eat9cROz0mwF0OZ/FsLIYaLroK9zNkLIYQQQ5wEeyGEEGKI67Vgr5T6k1KqTin1adKxPyulNsX/7FZKbYofr1BKtSed+23ScyYopTYrpbYrpR5SsiBcCCHEENCX0+i9mVX2JPAIsChxQGt9VeJrpdQvgJakx+/QWp/Vyev8BrgZWAcsAeYBS3u+uUIIIUTv0VqzefNmlixZwpIlS6ioqGDRokXHfmIP6LVgr7VepZSq6OxcvHf+L8Csrl5DKVUEuLTW78e/XwR8GQn2QgghBgGfz8fbb79tBPjkLcW/+OILotEoJlPvz6j31zr7C4GDWuttScdGKqU+BlqB/9RavwcUA8mbrVfHjwkhhBADjtaabdu2sXjxYpYsWcKqVasIBoPG+cLCQmMH09mzZ/dJoIf+C/bXAM8nfX8AKNNaNyilJgB/U0qdBnQ2P3/USQ6l1M3EhvwpKyvrweYKIYQQnevo6ODdd981eu87duwwzimlmDx5shHgzzrrrD4L8Mn6PNgrpSzAFcCExDGtdQAIxL/eqJTaAVQR68mXJD29BKg52mtrrR8HHofYOvseb7wQQggB7Nmzxwjub7/9Nu3t7ca57Oxs5s2bx4IFC5g7d+4RO432h/7o2c8G/qG1NobnlVJ5QKPWOqKUGgVUAju11o1KKa9S6nzgA+B64OF+aLMQQohhLBQKsWbNGmN4fsuWLSnnzz77bBYsWMAll1zCpEmTMJvN/dTSzvVasFdKPQ/MAHKVUtXAj7TWfwSuJnUIH2Aa8GOlVBiIAN/WWjfGz91KLLM/nVhiniTnCSGE6HUHDhxg6dKlLFmyhLfeeovW1lbjnNPpZM6cOSxYsIB58+YxYsSIfmzpsfVmNv41Rzl+QyfHXgJeOsrjNwCn92jjBqlVq1bx4IMPsnHjRmpqalJq3kMs6/O+++7jlVdeoaGhgbKyMr797W9z5513Go8JBALcddddPP/88ym18UtKSjp5RyFEb9BaEw6HMZlMA64HOJxFIhHWr19vDM9/9NFHKefHjRtnzL1PnToVm83WTy09frLr3QkKhUKsXLmS6upqSkpKmD59+lH3ue8pPp+P008/neuvv57rr7/+iPPf//73Wb58OU8//TQjR45k1apV/Ou//iu5ubl84xvfAOCOO+7g1Vdf5fnnnzd2vbv00kvZuHGj/NIRog9oramrq6Ompob09HTKysrIyMg4odeKRCK0tLQQCoVwOp0n/DrDWUNDA8uWLWPJkiW88cYbNDQ0GOfS09OZNWsWl1xyCfPnz6eioqL/GnqSJNifgM8//5zbbruNYDBIIBDAbrdjs9l49NFHGT16dK+9b+KOEkjp0SesXbuWb3zjG8bOeBUVFfzxj3/kgw8+4Bvf+AYtLS388Y9/5IknnuDiiy8G4Omnn6a8vJzly5czd+7cXmu7ECLG7/cba63b2tpoamo64SDd2NjI3r17AbBYLIwZMwa73d5jbR2KtNZ8/PHHRu/9gw8+IBqNGudHjRrFJZdcwoIFC5g+fTrp6en92NqeI8H+OAWDQW677Taam5uNY36/H7/fz2233cbSpUt7vYd/NBdccAH/93//x7e+9S1KS0tZu3YtmzZt4u677wZg48aNhEIh5syZYzyntLSUsWPHsnbtWgn2QvSDk6kAnpwBHg6HCYVCwzrYBwIBGhoaCIfDZGdn43A4AGhpaWH58uUsWbKEpUuXcuDAAeM5VquVWbNmGZ2pqqqqk/o3Gagk2B+nwwskJAsGg6xcuZLZs2f3catiHnroIb797W9TVlaGxRL7p3344Ye59NJLAaitrcVsNh+xDKSgoIDa2to+b68Qw1FGRgZlZWUcOHCA9PR0srKyTvi1nE4nhw4dAiAzM/OYgV5rTWNjIy0tLWRmZpKTk2P8rhgK6urqqKurQ2vNhg0b2LZtG2+88QarV68mHA4bjysuLjYy52fNmoXT6ezHVveNofOv3Eeqq6sJBAKdngsEAuzfv7+PW/RPDz/8MGvWrOG1116jvLycVatWcdddd1FRUcG8efOO+jyt9ZC8kxViIFJKkZeXR3Z2NiaT6aT+73k8HqqqqgiFQmRmZh5zVNHr9bJ7924AmpqaOr35H6z8fj9Lly7lrbfeYs2aNSm9d7PZzIUXXmj03s8444x+/Z3n9/tpbW3FYrHg8Xj65IZLgv1xKikpwW634/f7jzhnt9spLu6far7t7e3cd999/OUvf+Gyyy4DYPz48WzatIkHH3yQefPmUVhYSCQSob6+nry8POO5dXV1TJs2rV/aLcRw1RMJsUqp4+qVRiKRLr8fbHbs2GHMva9YsSKlI5aTk2P03ufMmXNSIyg9KRgMsnPnTqOtkUiEgoKCXn9fCfbHafr06dhstk6Dvc1mY/r06f3QqtjqgFAodMQvELPZbCSfTJgwAavVyltvvcW1114LxEYqtm7dypQpU/q8zUKIvpWZmYnL5aK1tRWbzTbohq/b2tpYvXo1y5YtY/HixXzxxRcp588991zmzJnDxRdfzNSpUwfkFEU4HE65KWlra+uT9x14V2KAs1qtPProo0fNxu/N5Dyfz8f27dsBiEaj7N27l02bNpGdnU1ZWRnTp0/n3nvvxeFwUF5ezsqVK1m0aBEPPPAAAG63m5tuuom7776b/Px8Y+nd+PHj+y3PQIiBwOfzUVdXh8lkIj8/f8guYbPZbIwcOZJAIIDVah0U68Srq6tZsmQJr776KitWrEhJSnS73cydO5dLLrmEuXPn9kkP+WTZbDY8Ho+R5O3xePrkfZXWQ7OE/MSJE/WGDRs6Pbd161bGjh17Uq+fWGe/f/9+iouL+2Sd/bvvvmssq0u2cOFCnnzySWpra7nvvvt48803aWxspLy8nG9961v8+7//uzE/1dHRwd13381zzz2XUlSntLS0V9veX3ri31oMbZFIhK1btxq9LYfDQWVlZb9sViJiPd/333/fGJ7/5JNPUs5XVlYya9YsrrvuOiZPnjwge+/HEgqF8Pv9mM1mMjMzeyx/QCm1UWs9sdNzEuzFUCb/1uJYgsEgn376KYnfhVarlXHjxnUaRPx+P+FwmIyMjEEZZAaquro63njjDZYsWcKyZctSljZnZmYya9Ysxo8fz5QpUygoKMDlclFZWdl/DR6gugr28mkVQgxrVquVkpIS9u3bB8SWZXUWyJubm42tSz0eD2VlZf1WU2Owi0ajbNy4kSVLlrB48WI2bNhAcsezqqrKSK678MILsdvtNDc3s3fvXiwWC0VFRf3Y+sFJgr0QYlhLLIVzOBwopUhLS+v0cfX19cbXzc3N5OfnS7A/Dk1NTbz55ptGYZtEfQCIrWSaOXMmCxYsYP78+Zx66qlHPN/j8eBwODCZTDLFcgIk2Ashhj2l1DGT8tLS0mhpaTEeL8P4XdNas3nzZmPufe3atSlL/crKyoyytDNnziQzM/OYrynX/MTJlRNCiG7Iy8tDKUUwGCQnJ2fI1EzvST6fj7ffftsI8Ik9ACAWqGfMmGEE+LFjx0oxrz4kwV4IIbqhP4tmJQuHw7S0tBCJRHC5XEeddoBYsa3m5maUUng8ni4feyK01nzxxRdGcD+8nHhhYaFRtW727Nm43e4efX/RfRLshRBiEKmvrzfKcqenp1NZWdlp7kAkEmHfvn14vV4gtpJg5MiRJ92bbm9vZ+XKlUaATyQtQmx64/zzz2f69OnMmDGD888/v8/WkYuuSbAXQohBJBG8IRZ4Q6FQp8E+Go2mVGfz+XxEo9ETKtO7Z88eI7i//fbbKYVtsrOzmTdvHgsWLGDu3LlEo1FjZcOuXbsYM2aMTHkMAL0W7JVSfwIuBeq01qfHj90P/CuQSMP8f1rrJfFz9wE3ARHge1rrZfHjE4AngXRgCXC7HqrFAYQYxNrb24lGo6SlpfVI3XfRObfbTWtrKxDb9e5oVfAsFgsFBQXGhjAFBQXd/ncJhUKsWbOGxYsXs2TJErZs2ZJy/pxzzjGG5ydNmpTyusnz9NFolFAoJMF+AOjNnv2TwCPAosOO/0pr/WDyAaXUOOBq4DRgBLBcKVWltY4AvwFuBtYRC/bzgKW92G4hxHFqaGgwdlPLz89nxIgREvB7SW5uLna7nXA4jNPpPGqGulKKwsJCo/79sbLdDxw4wNKlS1myZAlvvfWWcUMBsZuKOXPmsGDBAubNm8eIESOO+joul4tDhw4RjUZxu91DtvTwYNNrwV5rvUopVdHNh38JeEFrHQB2KaW2A5OUUrsBl9b6fQCl1CLgy0iwF2LACIVC7N271/i+rq4Oj8cz6DZZGSxMJlO3E91MJtNR/x0ikQjr1683huc/+uijlPPjxo0zeu9Tp07tdh19l8vFmDFjjB69LJcbGPqjMsF3lVKfKKX+pJRK7DlYDOxLekx1/Fhx/OvDjw87P/3pTzn33HNxuVzk5eVx2WWX8emnnx718TfffDNKKR58MGUQhUAgwL/927+Rm5tLZmYml19+ecqwmxBi6GpoaOC5557juuuuIz8/nylTpvDf//3ffPTRR6Snp3PJJZfw2GOPsWvXLj777DN+/vOfM3PmTMxmM+FwuNvvk56ejsvl6tWiQ4FAICXzX3Str2+5fgP8f4CO//0L4JtAZ+mhuovjnVJK3UxsyJ+ysrKTbetRaa1ZtWoVzz33HAcOHKCoqIhrr72WadOm9dq60XfffZfvfOc7nHvuuWit+eEPf8js2bPZsmUL2dnZKY/961//yocfftjpUNsdd9zBq6++yvPPP2/senfppZeyceNGGXYVJ8RqtVJeXs7u3bvRWlNQUCBDtwOE1pqPP/7Y6L1/8MEHxpbXAKNGjTLWvU+fPr3TufVEmVqlFBUVFQNixKaxsZHdu3ejlGLkyJGS8d8NfRrstdYHE18rpX4PvB7/thpI3natBKiJHy/p5PjRXv9x4HGIbYTTM60+4j34r//6r5SM1JqaGrZs2cLs2bP54Q9/2CsBf9myZSnfP/3007jdbtasWcNll11mHN+zZw+33347y5cvZ/78+SnPaWlp4Y9//CNPPPEEF198sfE65eXlLF++nLlz5/Z4u8XwkJ2dTUZGhpGgJ+VM+09LSwvLly83ytImEvQgdmM2a9YsY3i+qqqqy99XoVCI3bt3G5XvqqurGTNmTL8WwwmHw+zduxetNVpr9u3bh9PplM7KMfRpsFdKFWmtE5+8rwCJcejXgOeUUr8klqBXCazXWkeUUl6l1PnAB8D1wMN92ebDrVq16oilJxDLRF6+fDkzZ85k2rRpvd4Or9dLNBolKyvLOBYOh7nmmmv4z//8z053etu4cSOhUIg5c+YYx0pLSxk7dixr166VYC9OSk8XbBHdo7Vm69atxqYyq1evThlyLy4uNjaVmTVr1oDomZ8ss9ls3IBYLBapxNcNvbn07nlgBpCrlKoGfgTMUEqdRWwofjdwC4DW+jOl1IvAFiAM3BbPxAe4lX8uvVtKPyfnJfaB70x7ezvPPvtsnwT722+/nbPOOovJkycbx370ox+Rk5PDrbfe2ulzamtrMZvN5ObmphwvKCigtra2V9srRF8Kh8N0dHQQCATQWuNwOIbUzYjf7+edd94xhuf37NljnDObzVx44YVG7/2MM8444WBotVoZOXKkMYxfWlra74HVYrEwcuRI9u/fj1KKESNGyEhSN/RmNv41nRz+YxeP/wnwk06ObwBO78GmnZTkIbETOd8Tvv/977N69WpWr15tDF2tXLmSJ598kk2bNh3362mt+/0/sBA9JRgMsnv3bqxWK42NjQBkZGRw6qmnDupd6nbs2GEE9xUrVhAIBIxzeXl5zJ8/nwULFjBnzpyUEb+T5Xa7GTduHMCAGSp3OByMHj26v5sxqMiaiONUVFRETc1R0wZ6fZ/lO++8kxdeeIEVK1YwatQo4/iKFSuMZMGESCTCPffcw69//Wuqq6spLCwkEolQX19PXl6e8bi6uro+GY0Qoi/4fD7C4TChUMg45vf7CQQCgyrYBwIB3nvvPWN4/osvvkg5f+655xrD8xMmTOjV3u1ACfLixEmwP07XXnstW7Zs6XQoPz09neuuu67X3vv222/nhRde4N1332XMmDEp577zne/wta99LeXY3Llzueaaa/jXf/1XACZMmIDVauWtt97i2muvBWIJN1u3bmXKlCm91m4xvHR0dNDY2Eg4HMbtdvf55icmk4lgMIjT6aSjowMAm83W7XXi/am6utrovS9fvjyl3K3H42Hu3LlGWdqCgoJ+bKkYbCTYH6dp06Yxe/Zsli9fnhLw09PTmT17NhdeeGGvvO9tt93G008/zd/+9jeysrKMOXaHw4HD4SA/P5/8/PyU51itVgoLC43hLrfbzU033cTdd99Nfn6+sfRu/PjxzJ49u1faLYaXaDTK/v37aW5uBmKbtlRVVeFwOPqsDS6Xi6KiInw+HwUFBUYRmoEY7MPhMO+//74R4D/55JOU8+PHjzfm3idPniwFasQJk0/OcVJK8cMf/pCZM2fy7LPPGkPn1113HRdeeGGvzX0/9thjAFx00UUpx3/0ox9x//33d/t1fvWrX2GxWLjqqqtob2/noosuYtGiRTJMJ7oUiURobm4mEAiQmZmJy+Xq9LMeiURSyqxqrQkEAn0a7E0mEwUFBeTn5w/IXJS6ujreeOMNlixZwrJly4wbI4iVtJ09ezbz589n7ty5lJeXD8ifQQw+aqjuKTNx4kS9YcOGTs9t3bq106VpYuiRf+ueUV9fn5Lxfeqpp3Y6PK+1Zvfu3UZinFKK0aNHH7Mue7LW1lba2tqw2Wx4PJ4BeyMaDoeJRqNYrdYuA3I0GmXjxo3GpjIbNmwg+fduVVWVMfee6DDs37+flpYWsrOzGTFihPToRbcopTZqrSd2dk4+QUKIY/L5fCnft7e3dxrslVIUFxeTnp5OJBLB6XQeV6D3+Xxs377dWC6nlCIrK2vA9W79fj979uzB7/dTWFhIYWFhyk1JU1MTb775plHY5tChQ8Y5u93OzJkzWbBgAfPnz+fUU09Nee1Dhw4ZN0uHDh3C4XAcUSVTiOMlwV4IcUwZGRk0NDQY33e1Zt1ms1FYWHhC79PR0YHWmpycHLxeL7t37yYcDpOXlzegAn5zczN+vx+I1a9wOp3s3r3bmHtfu3atUfQFYuW7E2VpZ82a1WU54cNHW/t79FWW5g4NEuyFEMeUnZ2NUopAIEBGRkZKr769vZ329nbsdvtx9+I7OjpIS0sz5vTT0tKwWq0Eg0Fjk5N9+/bhcDgGVL19pRR+v5/169ezZs0aPvjgg5QluRaLhRkzZhgBfuzYsd0OmIn96hPD+C6Xq7d+jC6FQiFqa2tpbW0lLy+P3NxcKV4ziEmwF0Ick8ViSanNkNDW1sa2bduIRCIopTjllFO6tdSutbXVGK5XSlFZWYnT6cThcDBq1KgBWdFRa80XX3zBkiVLeP3113nvvfdS1vIXFhYamfOzZ88+4SWHdrudkSNHEolEsFgs/RZgm5ubqaurA2I3XGlpaf124yFOngR7IUSXgsEgbW1tKKWO2HDE5/MZw9Vaa7xeb7eCXFtbmzE8rbXG7/cbNdsdDgdFRUVG7z6RA9Af2tvbWblypTE8v2PHDuOcUorzzz/f6L2fddZZPRaYzWZzvycmJu+O19n3YnCRYC+EOKpgMMiuXbuMBL3CwkJGjBhhDEkfXpGuu1njh695P/z7zMxMRo8ebWS796Xdu3ezdOlSFi9ezDvvvJNSTyM7O5t58+YZhW0O32eiu5LnwbXWhEIhTCbTgMq6d7vdNDc34/P5yM3NNaZaotGokUTodrvJy8sbVJUJh6uB88kSQgw47e3tKZn4tbW15ObmYrfbgVhAKC4uprGxEafT2e2scY/HQ3l5Oe3t7UfkACT0Ve82FAqxevVqo/e+ZcuWlPNnn3220XufNGnSSbWpvb2dAwcO0NHRQVFRER6Ph0OHDlFdXY3FYmHUqFF9WpOgK2lpaZx66qlHTCd4vV6qq6uB2KoEu91OTk5OfzZVdIMEeyHEUR0e2KxWa8oxs9lsLD073tc90V7x0XR0dNDR0YHdbj/msP+BAwdYunQpS5Ys4c0338Tr9RrnnE4ns2fPZtKkSZx11lmUlZUxcuTIHplKOHToEE1NTQDs2rWLyspK9u3bB8RuOurq6gZMsIfOb7gOH85PXnUgBi4J9kKIo8rMzKS8vJz9+/djNpspLy/vtaHmSCRywr1mv9/Ptm3bCIfDmM1mKisrU1YGRCIR1q9fb2wq8/HHH6c8f9y4cUZy3dSpU2lra2Pnzp1ArDfe2traI8E+OTBqrY3h/ET+wmDIdnc4HGRlZdHU1ITD4ZCkvUFCgr0Q4qiUUuTm5hpbpp5oMPb5fLS0tGA2m8nKyjKmASDWoz148CBNTU1GXfvjrWPf1tZGOBwGYgHV6/XS0dHBsmXLWLx4MW+88YZRqAZie1nMmjWLSy65hPnz51NRUZHyeokNdBJ6ajohLy+P1tZWwuEwxcXFOJ1OKioqqKmpwW63H7G/xUBktVopLy+nuLgYs9k8oPIMxNHJv9JJ8Pl8NDc34/F4+mTo7dFHH+V3v/sdu3fvBuC0007jP//zP7nkkkuAWE/hv/7rv3j88cdpamrivPPO49FHH+W0004zXiMQCHDXXXfx/PPPG7XxH3vsMUpKSnq9/WLgiEajx9WLPJlg19HRwfbt241ebXt7OxUVFUaCWmtrKwcPHgRiZXnT0tKOe0c3i8VCNBrliy++YPXq1WzcuJGNGzemDDmPGjXKmHufPn16lz11h8NBWVmZkYvQUzv3ORwOxo0bRyQSwW63o5QiOzsbj8eDUqpfitecSNGcgbBaQBwfCfYnYP/+/fzyl79k7dq1WCwWwuEwU6ZM4fvf/z7FxcW99r4lJSX87Gc/o7Kykmg0ylNPPcWXv/xlNm7cyPjx43nggQf4xS9+wZNPPsno0aP58Y9/zMUXX8znn39uLGu64447ePXVV3n++eeNXe8uvfRSNm7cKP95h4FgMEhdXR1NTU14PB7y8/NTetm99Z7Jw9der9dI+gKMHnnC4d8fTTgcprW1lRUrVrB48WIWL15srAuHWA901qxZRt35ysrKbgc1k8lEXl6eUVvA6/VSX1/f7Xr90WjUyANwOBwpj7darUdkr/fH8L3WmoaGBg4ePGgsdxyIOwOKntFrG+Eopf4EXArUaa1Pjx/7OXAZEAR2ADdqrZuVUhXAVuDz+NPXaa2/HX/OBOBJIB1YAtyuu9Ho3toIZ//+/Xz961+nra0tpddgMpnIzMzkmWee6dWAf7js7Gx++tOfcvPNNzNixAi++93v8oMf/ACI9aDy8/N58MEHueWWW2hpaSEvL48nnniC6667DogVyygvL2fp0qXMnTu3z9rdV2QjnFQHDhxIqfRWUFDQo6M6fr8fv9+P1Wo1dsYLBoNs27bNGBo//D3b29vZvn07gFFYx+PxpAwPh0IhotEoFouFTZs28dprr7FkyRI2bdqUcnNQXFxsBPdZs2YZN7kn+zN9/vnnxv/38vLyYyYXJl/nw5crDhQ+n4/PP//c+L64uPiEyxyLgaG/NsJ5EngEWJR07C3gPq11WCn1M+A+4J74uR1a67M6eZ3fADcD64gF+3nA0l5q8zH98pe/PCLQQ+xOvq2tjV/96lc8+OCDvd6OSCTCX/7yF3w+H1OmTGHXrl3U1tYyZ84c4zHp6elMmzaNtWvXcsstt7Bx40ZCoVDKY0pLSxk7dixr164dksFepDp8Q5tE79hms6XUuz+Rod3kanrwz6Bos9kYNWoUPp/P2Fs+WXp6OpWVlTQ0NFBbW0tDQwOtra2UlZVhsVioq6vjxRdfZOXKlaxbt85Y9gWx4eRzzjmHr371q1x66aWcccYZPR5UEzcaCYFAoMvHh8PhlI1vDh06RH5+/oBbi354n0mK5gxtvRbstdar4j325GNvJn27DvhaV6+hlCoCXFrr9+PfLwK+TD8Fe5/Px9q1a4/6nyIajbJmzRp8Pl+vzeFv3ryZyZMn09HRgcPh4JVXXuGMM85g7dq1AEfMdRYUFLB//34gtka6syVPBQUFA7I8qeh5ibrrEBsVam1tpaGhAaUUo0aNwmazUVNTQ0dHB/n5+cdVD72trS1luL6+vt74rKWnp3c5R242m1MC5ObNm3n55Zd5++23eeedd4w6+Yl2T548malTp3L++eeTnZ3Naaed1muJYom2t7e3o5Q65v9ts9mM0+k0EgIPH8YfKDIyMigoKODgwYNkZGQYSZhiaOrPOftvAn9O+n6kUupjoBX4T631e0AxUJ30mOr4sX7R3NyMxWJJqYd9OIvFQnNzc68F+9GjR7Np0yaam5t56aWXWLhwIe+++65x/vBeTXd6aLKr1fCR2NDG5/NhNpuNz3Ji/jYcDhu9/3379mG327udnHZ4z7WrnfEOFwqF+Oijj1i+fDlr1qxh7969KefHjRvHBRdcwIIFC6ioqMDpdBrr1UtKSno1I9xms3HKKafQ3t6O1Wo95mY/iW1+MzMz0Vrj8XgG5JI6s9lMcXExeXl5klU/DPTLv65S6gdAGHg2fugAUKa1bojP0f9NKXUa0FkEOup8vVLqZmJD/pSVlfVso4lV/TpW8lA4HMbj8fT4eyfYbDZj/+uJEyfy4Ycf8qtf/cqYp6+traW0tNR4fF1dndHbLywsJBKJUF9fn7KpSV1dHdOmTeu1NouBI7GhTV5eXspwOMTyTg5fcnb4ja3X66Wurs5IYEu+qXW5XBQXF1NfX09GRsYx53/37dtnFLZZvnw5bW1txjm3222UpZ06daqR1FdYWEg0GqWjo4Py8nIyMjJ6dTe8RKJdNBolMzOz2wlsNpttUCyjU0r1eoKmGBi6FeyVUq8BLwCvaq3bjvX4Y7zWQmKJexclEu201gEgEP96o1JqB1BFrCefnD1UAtRwFFrrx4HHIZagdzLt7IzD4WDKlCm89957nQ7lm0wmpk6d2qcVsKLRKIFAgJEjR1JYWMhbb73FueeeC8SWPL333nv8/Oc/B2DChAlYrVbeeustrr32WgCqq6vZunUrU6ZM6bM2i4EhUao18Vl2u93YbDZjSsdsNqcE0kQiXeLxra2tjBkzxggWiWp6+fn5nfZkw+Ew77//vlGW9pNPPkk5P378eKOwzeTJk1N6momMfpvN1qdD4vX19UaFO7fbTUVFhfSAxaDU3U/tL4CrgJ8qpdYTG35/XWvd0fXTUiml5hFLyJuutfYnHc8DGrXWEaXUKKAS2Km1blRKeZVS5wMfANcDDx/Pe/a073//+3z00UedZuM7HA7uvPPOXnvve++9l0suuYTS0lK8Xi/PPfcc7777LosXL0YpxR133MFPfvITxowZQ1VVFf/93/+Nw+EwArvb7eamm27i7rvvJj8/31h6N378eGbPnt1r7Rbdl6hFn5j3PXxoPBqNEgwGjyhbeyIcDgdjx47F7/djs9nIzMzE5XKRnp5OOBw+Yg/5wxPVwuEw4XD4iJ5hcqCvq6vjjTfeYMmSJSxbtozm5mbjXGZmJrNnzzYCfFerAjrrUYfDYRoaGvD5fLjdbrKzs09ouDwajRIOhzvdTrahocH4uqWlhVAoJMFeDErd+tRqrVcCK5VSZmAW8K/An4Cj1klUSj0PzABylVLVwI+IZd/bgbfic8SJJXbTgB8rpcJABPi21jpR7upW/rn0bin9mIkPseUpzzzzzBHr7KdOncqdd97Zq8vuamtr+frXv05tbS1ut5vx48enLJn7j//4D9rb27ntttuMojpvvvlmyvKjX/3qV1gsFq666iqjqM6iRYsGZALRcHN4NrvL5WLkyJFGcPH7/VRXV+P1esnIyKC0tPSkR5HS0tJS5tYtFstRN7Ox2+1YrVZjaD89Pb3Tm5ENGzYYvfcNGzakZH2PHj2aefPmcemll3LhhRee1BByS0sL1dXVmM1mlFJYLBbcbvdx5Z90dHSwb98+vF4vOTk5FBcXpwRzl8uF3+83fl4J9GKw6vY6e6VUOrE18lcB5xDr2f9bL7btpPTWOvtkfV1BTxy/wbTO/vA18ACnnnoqbrcbrTW7d+9OKfmamZlJZWVlpzdqiS1Tj+cmrr293Vgj73Q6Ow2abW1ttLS0oJTC4/GQnp5OU1MTb775JkuWLGHp0qUpWfV2u50JEyYwZcoULrvsMoqLi2lrayM7O5vCwsJuBc9IJILJZDqiPQcPHuTAgQNkZWVRX18PxHJ1kvNRjuXgwYMpuQuJXefMZjMmk4lQKERzczPRaNQY9RBioDrpdfZKqT8D5wFvEFs7v1JrPewXZTocDgnyosd0deOdqPeeYDabSUtLo76+HrvdjsPhMEaZDh06xMGDB7HZbBQXF3crm/7wUYXS0tJOE8wyMzPJyMhg8+bNPPXUUyxevJj3338/ZcldWVmZUZa2tLSUcDhs3BwkMugPHjxIenp6l1ujaq2NnyUtLY0RI0YYm+VkZGTgdDppb283Aj3Ekv5cLle3RwySr3lmZiZNTU3s3r0bj8dDcXExNpvtuG4ehBioujsm9RSwHrgA+C5wplLqt8c7Zy+EODqXy0Vtba0RgBKBFWLD61lZWUY52Ozs7JQedF5eHqWlpTQ3NxujA+3t7ezcuZOxY8cesQwuHA7j9XrRWuNwOIxs94Samhqys7MJBAKEw2FCoRBr1qxh8eLFLFmyxKjdkGjbjBkzjAA/duxYoxe+Y8cOmpubMZlMx10W1+fzGclxZrOZmpoao0bAyJEjyc7OJj8/n6amJiOXwGQyHde8vcfjwev14vV6yczMNK5voiZ+T23DGw6H0VoPuMI6YvjobrC/gdj694fi318DPA1c2QttEmJYcjgcjB49+qgJevn5+SiljC1Xkx06dIjs7Gza29tTjidWayQH+3A4zL59+1KKvhze+zebzWzYsIEXX3yRNWvW8PHHH6cswyssLDQS62bPnn3U0YPCwkJjn3mn02mMTpjN5mOOiiUnA6alpRmjAhAbGcjKyiIjI4NRo0axZ88eIFa1r6uA6vf78Xq9xvx+Wloao0aNIhKJHHFNe4rX62XXrl1EIhHKy8uPmhMhRG/qbrAfrbU+M+n7FUqpv/dGg4QYDhJZ9RaLJWXeOjMz86hFW+x2OyUlJQSDQbZu3ZpyLlHMKbl3DrGe7uFD2h0dHSlz/z6fj4KCAqxWKytXrmTt2rWsW7fO2F0RYuuxJ06cyJe+9CUWLFjAWWed1a0edGZmJqNHjzY2vnE6nQSDQdLS0sjIyMDn83Hw4EG01hQWFqbcAGRmZpKTk2NU+Eve9z09Pd0YPXC73Zx22mlorbvMAQgEAuzYscOoxpeoBZ/YwS2R0Z/Iw+mpfdpramqMG6Xdu3eTmZl53ImJwWDQyB3weDzHVbBoKNFa09zcTHt7O5mZmT22G+Fw0N1g/7FS6nyt9ToApdR5wJrea5YQQ5ff76empoaWlhasVitlZWXHVYjJZrOllEEGyMnJSZkHb2lpMebsDw8MyUG6pqaGNWvWsGnTJlatWpUyMuB2u5k8eTIXXHAB559/PmecccYJFYpJvqFJDuahUIidO3cagbCtrY1x48YZPXOLxWLkDpjNZmPqwm63HzGP3p1ExGAwmFJ2t7W1NaXwT2Kf9pKSEiNBrzOtra00NjZit9vJyck5ZqGd5Nc50W1sa2trjWmblpYWTj311GG5gqalpYWdO3ca348ePVryprqpu8H+POB6pVSihmUZsFUptRnQWuvxvdI6IYYYrTUHDhygpaUF+GfAGzt27HFleufm5mKxWGhrazOWzx08eJD29nZj3wSXy3VE7zQUCrFu3Tpefvllli5dmvKLE+Ccc85h5syZnHbaacyYMcPoTSbW4PekSCSSMjUQDoeJRCIpw/DJhX26W7o3sVGNzWZLCaw2mw273W6c7+y1jjXn39HRwY4dO4wpBqXUMSsFJhILw+EwpaWlx72NrNY6JTnT5/MRDodPKthrrfH5fMa/7WBZUnj4JkTJN2+ia939F57Xq60QYphILOVKprWmvb39uIK9xWIhNzeX3Nxc/H4/u3btwu1209LSYvwiLyoqAmJL+hJlad98882UwOF0OpkzZw4LFixg3rx5jBgxgkOHDrF3714aGxvxeDxGQDveoeOmpibq6+tJS0sjLy/viOcnSsomkuISO+SdjObmZnbu3InWmvLycnJycoyAb7fbOeWUU1Lm7I9XJBJJySVIDjaBQIDW1laUUrjdbuOmJTGVobU+6o1ENBqlubmZYDBIZmZmSm0MpRR5eXlGsmJP7KDX2NhoTNPk5ORQWlo6KEYKEssiI5EIdru9V0slDzXdLaqzp7cbIsRwkOipJgq1JBw+h5vo5R7eO+1MS0sLHR0dKKXIzs4mHA6zZ88eHnjgARYvXszHH3+c8vhx48YZe75PmTLliACbmZlp/EJtbm4mPz//uAN9W1ubMWrQ2tqK1vqI/SpMJhMjRowwRgwyMzNPasOYaDRKdXW1Ma+/Z88enE5nyrU91u57x5KWlmbcoJjNZmOnuEgkwt69e40kv7y8vJSf91jD983NzezatQuIXZfRo0enBLK8vDzS09ONnvjJbqyTvJKjoaGBgoKCQVFDIHHjlMj7kLr+3Tc4xm6EGCISO43t3LnTSKYrLi5O+cXu9XrZt28f7e3tuN1uiouLU34RJwrbJAq9JJKWli5dypo1a3j//feNaQKIBbiLLrqIBQsWMH/+fCoqKrpsY0ZGBlVVVbS1tWE2m3G5XMcdXA5fVnf4BjvJ16Mnk6ySe6edFeLpidcfMWIEOTk5mM1mI9iEw+GUbP6WlhajJkB3JF+fxAqK5M+EUiqlt3+yHA6HsfGQzWYbNMP4cPI3bMPV4PkXFkPKjBkzOP3003nkkUf6uyl9InkduMvlYuzYscaWqRkZGUZQCgaD7Ny50wiWLS0tmM1mKioqUEoRCATYvn07wWCQL774go8//pgVK1bw8ccfpxSIGTVqlLHuffr06V3+cgyFQkZwcjqdxk5yJzNEmp6ejsPhMLbLTU6oCwaD+Hw+I4D1VKAxmUyUlpayb98+IpHICc2Pd8fhGwTBP+sgJJYHZmVlHdewePIKDKvV2uvBLDEVEIlE8Hg8sv5/GJBgL/rFyy+/PCh+wSS2U9Vak5aWdtzzmoFAgJaWFurr64lGoxQUFJCdnY3dbsdut6O1prW1FZ/Ph8ViwW63H9ErbmlpIRwO4/f7+dvf/sbLL7/MmjVrUjZpsVqtTJo0iblz53LllVcyevToI3q1fr/f6H263W5jWHj//v3Ga1ksFkaPHn3SS7tsNhsjR47E7/djsViMYJaYYki0Iz8/n+Li4h7b793hcDBmzJgu58d7g9lspqSkxKjNf7zJjG632xieTk9P7/WldYkVHWL4kGA/iAWDwV7pufTW6yYbDIVFOjo6qKmpMXprTqeTkpKSLnu8fr/f2IktJyeHjo6OlPnRvXv3YjKZjDKxiYSyhLy8PGO+XGvNzp07+eijj/jggw9YvXp1yo1AYWEhkydPZurUqZx77rlkZmZSWFjY6WZMgUCAbdu2Gc8/dOgQVVVVQOrObuFwmPb29h4JNjab7YjPUSKJLaGuro78/PwenXs90eVtJ8tms3VZ/vdYZAmZ6E19d+srTtqMGTO49dZbueuuu8jLy2Pq1Kls2bKFSy65BKfTSX5+Ptdcc42xHznEfnnfeeedZGVlkZWVxZ133smtt97KjBkzunxd4JivvXnzZi666CJcLhdOp5MzzzyTFStWALGh4e9973uMGDECu91OaWkp9957b8p7fve73zW+b2pqYuHChWRlZZGens7s2bP57LPPjPNPPvkkDoeDt99+m9NPP53MzExmzpxpJDX1NK01Bw8eTKna5vV6qaurM/Y4r62tTUm0CwaD7Nixg7q6OiwWC7W1tSlLy5J/1sR7JN8IQKy2+5YtW/jZz37Gl770Ja666ip+9rOf8e6776K15sILL+S+++7j+eefZ9WqVfzwhz9kxowZRs/58BsRv99PXV0dTU1NKTcKwWAQv9+P2Ww+Ysi4N0dcLBZLyujIiYyW9LZwOExjYyMNDQ2d/vsJMRhJzx76pRcAXW98cjTPPPMMN998M++99x5NTU1MmzaNm266iQcffJBQKMQPfvADLr/8ctatW4fJZOLBBx/kySef5A9/+AOnn346jz32GM899xxnn332UV83sRb8WK997bXXcuaZZ7J+/XosFgubN282eoQPPfQQr7zyCi+88AIVFRVUV1fz+eefH/XnuuGGG/j888959dVXycrK4gc/+AHz5s3jiy++MIJRIBDgpz/9KX/6059IS0tj4cKFfPvb32bZsmXHfR2PJRgMpmywArGeW2IYOuHgwYNUVlaSkZFBW1ubsRTLbDYTCoU6/WxZrVbjZgKgurqa1atXs2bNGjZu3JiynCsrK4spU6Zw6aWX8rWvfY3c3FxjqZ7WmkgkQl1dHYFAgIKCgpRkt/b2dr744gsikUinIymJYjfl5eUcPHiQUChEfn7+MXuYiaV9HR0dpKWl4XA4uv1/KLH8LZHNnpeXN6CSwxKf/cRywKysLCoqKo45JeDz+Ywh+OQSxbI0TAwUA+d/meiWkSNH8otf/AKAH/7wh5x55pn87Gc/M84vWrSI7OxsNmzYwKRJk/jf//1f7rnnHr761a8C8Otf/7rT4Jj8ut197T179nDXXXcxZswYILYda8KePXuoqqriwgsvRClFWVkZU6ZM6fRn2rZtG6+99horV65k2rRpADz99NOUlZXx7LPP8q1vfQuI9bgeffRRRo8eDcBdd93FjTfeSDQa7fKXcSI4tba2EolEjDKbXQUZpRQmk8lIrDObzXg8HiMIJCRvKJPcC4xGo0YVtuTXSeyJ/sorr/DXv/6VDz744IjCNhMnTmTmzJmcfvrpjB071kjqC4fD1NfX43A4UoJIIkv78IDr8/mMjH+fz0d2drYxqjBixAhjNCAzM5NRo0Yd9VocrqWlhR07dgCxnnlpaelxzVE7nc4ezSzvSdFoNKWUcFNTEyUlJV1Oa7W2trJ9+3a01uTm5tLY2Gh8Jquqqo5a/liIviTBnhPrYfeXCRMmGF9v3LiRVatWddoT27FjB6NHj6a2tpZJkyYZx5VSnHvuuUaBjs5etzuvPWnSJL7//e/zrW99i6eeeoqLLrqIr371q0bgv+GGG7j44oupqqoyirbMnz+/06C8detWTCYTkydPNo653W7OOOMMtmzZYhyz2+1GoIdYwEoUqekqB6CpqSlluP/QoUNkZWVRVlZGJBLB7/cTDodJS0sz1jAnCr7U1taSnZ1NJBI5ou48xIJde3s71dXVOJ1ObDYbwWCQ1tZW4xe/x+OhpqaGtWvXsmzZMtauXZtSltblcnHhhRdy8cUX89WvfpWSkhKqq6uNnr/T6UyZ605LS6OystIIQEfrVScPjweDQaPMqs1mO6k5+eQM+8TOdInktMHOZDLhdruNPAaXy3XMaYa2tjbjd0hy0Z1oNIrf75dgLwYECfaDTPIvjmg0yiWXXMKDDz54xOMKCgpSSnoez+t257UB7r//fq677jqWLl3KsmXL+K//+i9++9vf8s1vfpNzzjmH3bt388Ybb/DOO++wcOFCzjzzTN56660jAn5XN1vJbT+8J544l1zR7HBa6yNubDweD3a7nUOHDtHY2JiyxrmgoICioqKUIeba2lrC4TAejyellw6x65YIDF6vF4/HQ2ZmJsFgkI8//pj33nuPZcuW8emnn6a0YfTo0Sk15y0WC6NGjTKKtCSmLpRSZGRkpFS96+joMHrqXXE6nUZVPbvdTkFBgVEw52RYrVbcbjeNjY3GDdCOHTsYN27coN+gRSlljHporXG73ce8Xsk/8+Gf7d5OdBWiuyTYD2LnnHMOL774YpfbehYWFrJ+/XpmzpwJxILfhx9+eMx63t15bYDKykoqKyv53ve+x6233sof/vAHvvnNbwKxYHPllVdy5ZVXcsMNN3D++eezfft2Iws8Ydy4cUSjUd5//31jGL+1tZXNmzdz4403dvt6dEZrnZKYlpWVRSgUwufzYbVa6ejowGQy4fF4jKH4RNC22WyYzeaUNe+J3noiyCVuOEwmE5FIhJdffpmVK1eyevXqlMI2mZmZXHDBBUyYMIGpU6dSUFBAVlaWMV3gcrlSNsNxu93GyEVn+753Z1mZ1Wpl5MiReL1eDh48yN69e/F6vRQVFZ3UOu7EDUnyz5fIHxgKbDbbERvtdMXj8XDKKacYhXDcbrdR/rin9xMQ4kT1WrBXSv0JuBSo01qfHj+WDfwZqAB2A/+itW6Kn7sPuAmIAN/TWi+LH58APAmkA0uA2/VgGnfvRbfddhu///3vueqqq7jnnnvIy8tj586dvPjii/ziF7/A6XRy++2388ADD1BVVcW4ceP43e9+x4EDB4y66Sf62haLhbvuuosrr7ySiooKDh48yOrVqznvvPMA+OUvf0lRURFnnXUWVquV5557DpfLRUlJyRHvVVlZyZe+9CVuueUWHn/8cTweDz/4wQ9wuVxce+21J3WNkrdGNZlMaK2x2Wy0t7cbc+w5OTkpWfGJ4W6Hw5ESwLTW1NfX43K5jF/ka9eu5aWXXmL9+vX8/e+puz5XVVUZhW0uvPBC/H5/yvx8U1MTRUVFjBgx4oh2WywW8vLy+Mc//oHJZErZCz4rK6vbQ8NKKerr642h96amJiwWyxGla49HIhgGg8GURLaBWro0EAjQ1taGzWbrleVtSqkjdi1M3BAJMVD0Zs/+SeARYFHSsXuBt7XW/6OUujf+/T1KqXHA1cBpwAhguVKqSmsdAX4D3AysIxbs5wFLe7Hdg8aIESNYs2YN9913H/PmzaOjo4OysjLmzJlj/OK96667qK2t5cYbb0QpxY033shXvvIVYz74ZF47sVyutraWnJwcLr30UmPY3+l08vOf/5xt27ahlOLss89m6dKlR81OfuKJJ7jjjju4/PLL6ejoYOrUqbzxxhvH1QNNvgfUWhvf5+XlUVdXh91up6Ojg/T0dKPuvN1uT1leB7F514aGhk4Ls7S0tLBs2TI2bdrEihUrUm4SbDYbEydOZMqUKUydOpXJkyenFC6xWCyMGDGCAwcOoLUmKyury3XZkUjE2OUrLS2N7OxslFIUFRV1e3lcJBJJWdcOGAmFJ7MKJbmuvdYah8MxoLLqExIVCRNLJE855ZTj2k54MEqMXCVuEvuyuJAYuFRvdpKVUhXA60k9+8+BGVrrA0qpIuBdrfXoeK8erfVP449bBtxPrPe/Qms9Jn78mvjzbznWe0+cOFFv2LCh03Nbt25l7NixJ/nTDV7nnHMOU6dO5eGHH+7vppyUaDRqFJ9RShlfJ4bUtdZs377dyF9oaGggPT0dk8lEc3MzOTk5RCIRI9gnlkolito0NjZisVioqalhyZIlrF69ms2bN6f09svKypg3bx7jx4/nrLPOIi0tDZPJRFZWFiaTCavVitPpTOlRBgIBotEoaWlpXQZcrTW7d+9OyQ7Pz8+ntLS029dIa82ePXtSCuccbTRhKGptbWXbtm3G9wUFBZ2OLg0VkUiEPXv2GJ/p0tJS8vPz+7lVoq8opTZqrSd2dq6vb8ULtNYHAOIBP/EpLCbWc0+ojh8Lxb8+/HinlFI3ExsFOKlhyqFkz549LFu2jOnTpxMOh3n88cf5+9//zuOPP97fTTspiV6v1hqLxWIMyVsslpQ5bqWUsSa9sLDQWD+fk5NDc3MzBQUFtLS0GDcBjY2NZGZmsnjxYtasWcOaNWtSlttZLBZmzpzJggULWLBgASNHjjT2kD9w4AAQqw6YvEZfKUVlZaWx3Oxow90dHR3GqgCLxWL04u12O16vF5fLddwV2hLb01osFlpbW485mjDUWK1W4+YNGPQJhMcSCoVSRqoaGxsl2Atg4CTodda90V0c75TW+nHgcYj17HumaYObyWRi0aJF3H333USjUcaNG8fSpUuZOLHTm79BIZF0l+jFJ28yk5wpbzab0VrT2NhIVlYW1dXVRKNRsrOz0VqTn59v1KvfunUrixcv5sMPP+SDDz5IWTOfk5PD1KlTueyyy7jyyiuNJWatra384x//IBqNkpGRQV5eHm1tbcbQe3J7m5ubu1xbnqjKl3it8vJyMjIySEtLO+leeFpa2pDuzXYlPT2dU089Fa/Xi81mw+PxdLl3/GBnsVhSNiCSBEGR0NfB/qBSqihpGD/RZaoGkscmS4Ca+PGSTo6LbiotLWX16tX93YwelZxhbzKZjK+VUinBPnmIPPlcYrndRx99xPvvv8+6detS1uErpRg/fjxTpkzhggsuoKqqCpPJREVFhRHoI5EINTU1xmv6/X78fj+VlZXs37//iDYnphXC4bBRbCehvb2dvXv3GjkGfr+f+vp6GZ3qIQ6HI2Uapb6+3qiC2Nne8YNZoiJia2trj28fLAa3vg72rwELgf+J//1q0vHnlFK/JJagVwms11pHlFJepdT5wAfA9cDgnmgWXTp8Ht5kMh2xzjk5wz7xOK21UbUsEokYxxK01tTU1BhD8x9++GFKD9zj8TBt2jQmTZrERRddhNPpTFl/n56entIDTBTjOZzf7yc3N5e9e/emHHc6nezevZvm5mYcDgeFhYXG64VCoSNqDfREEp3o3LH2jh/s0tLShvx0hTh+vbn07nlgBpCrlKoGfkQsyL+olLoJ2AtcCaC1/kwp9SKwBQgDt8Uz8QFu5Z9L75bSQ5n48ot04EmehwdS5lqTA75SCqvVSjAYJBKJGPP0ycP6iWV2gUCARx55hHXr1qUkagGMGTOGWbNmce2113LGGWdQX19Pc3OzsTwvIyPD+Jzk5eWlFEixWq14PJ4jMvkTteKj0SgHDx5EKUVxcTEHDx40bg5aW1tpa2tj7Nix2O12Y44+OdcgkXl/ooLBIE1NTQQCARwOh1EQSKQWkLJYLCcVGBOVDS0WS7eq7QnRX3ot2GutrznKqYuO8vifAD/p5PgG4PQebJqxSYlUtxo4EgVttNaYzWYjux7+WYI00aNPlK61WCzGucSyr/b2dtrb22ltbcXr9bJ7926efvppIPZL/rzzzmPq1KlMmTKFwsJCTjnllJTa8iaTiYaGhpTlaomKai0tLbS0tBCNRklPTzd2jGtvbzcqr7lcLkwmEwUFBeTk5BgJgoePAiRGBux2OzabjVGjRlFbW0t7ezs5OTknlUSXGMVIZOAfOnSIioqKYZWY1xWPx8Opp55KKBQiIyPjhAsMhcNhdu/ebcyPl5SU9Pse8Yn/G1arVTozIsVASdDrUx6Ph4MHD1JcXCy9nQEiMXyf6JEnAn2ix6uUMhLuEj3gxGYzfr/f2OgmeYi2ra2NNWvWcOONN3LeeecxZcoU8vLyjJsDp9OZ8os+8X1GRgaHDh0yevRZWVkp9fVzcnKoro4tEklPTzcKyuTn56d8nhI3IEfr7SUfTyzPS/xMJyMUCqUstYPYaIIE+xilVI/MZYfDYSPQQ+wa92ew7+joYM+ePbS1tVFQUEBhYaGMNAjDsAz2ubm5x9xyVfStaDRKOBxOyahPzLsnz9ErpQiHwwQCAQKBAB0dHSnz3Uop7HY7VqsVq9XK9ddfT0ZGBna73ZjL7yq73WKxkJ+fT25uLvDPBMCamhrj9ZPX2SdGEgCjLv7h0tLSKCwspLa21jjmdruPmCdO3NCcLJPJRFpaWsqNj8zh9rzE0H1iFKi/k+Gam5uNm4/a2lqcTqdk4wvDsAz2JpNJMp37SKJUaaKIjM1mo6Ojwxj6djgcpKenGwVwLBaLMQ+enZ1t7Bq3ZcsWlixZwgcffMCmTZtSAnxxcTEXXXQREydO5JxzzqGoqChlLj2RKQ90u+Z5cg89cXNxohJr3TMyMujo6MBms+FyuXqt4lyiHO6ePXsIBAK43W4p39oLEpnvXq8Xs9nc74H18GF7GcYXyYZlsBc9x+/34/V6CQaDpKenGwlqECunu2fPHiKRiFFLvqOjI2UDFZvNRnFxMRkZGRQVFRnFaHw+H+vWreOtt95i7dq1KWVpLRaLUQXwggsu4LTTTsPpdBpFbRLFYw5PnvN4PCdUG91qtZKRkYHf7zdyCg7n8Xi6nPs1m819GnCdTidjxowhHA5js9mG/XSV1pq2tjZCoRCZmZk9lq9js9kGzPRIVlYW7e3t+Hw+8vPze2UfADF4SbAXJ6ylpYUdO3YYw5eJevJZWVm43W5j05fEsHo4HE4J9Dk5OQQCAXbt2oXWmtbWVlatWsXrr7/O3//+95Ts9MLCQqZMmcKMGTO44IILjIS93Nxco0BKIsBHIhEj0S2Rpe/xePB4PCfUmzabzYwYMYIdO3agtTZ2v/P7/QSDQXJycsjNzR1wAdVisQzIevX9obm52fg8OhwORo4cOeQSdG02GxUVFcZqFCGSyW8CcUKCwSC7d+/G6XTS3t6eMszd2Nh4xJ7vbW1tRywtq6mpYf369axZs4Z169YZSW8QC7ATJ07kggsuYM6cOcY8en19vbEszm63U1tbawzpK6WMPd5tNptRHa0n5sHdbjdVVVX4fD4ikQiZmZmMGDECpZQE1EEgeZTH5/MZ0ylDkQR60Rn5LSVOiN/vJxwOGzXXk9lstiP2YE8Uu6murmbDhg2sWrWKdevWEQwGjcdkZWUxZcoULrroIs4991zS09NT1oe3tLQYvfVAIIDVak2Zu08k7QGMGjWqx3c3O7wSmxg80tPTjYCf2KBIiOFEgr04bsFgkEAgYATwZFlZWaSnpxsJcVarldWrV/PGG2/w/vvvp5SlBTjttNOYPXs255xzDmPHjjV2jEv8Ym5ubjbKfubm5hIKhYwh+eR5/MMd3i4xvCWmWcLhMC6Xy8ivSCz3lGQ2MdRJsBfHFI1GjXXsZrOZQ4cO4XQ6j6hCmJWVRTQa5bPPPmP9+vW8+eabfPjhhykFZVwuF1OmTGH69Omcc845nHbaaTQ3N6csE2tvbzcS4gBju9mEvLw8LBZLlwlxJ1ooRQxNVqv1iDXwzc3N7Nu3D6vVSmlpaafLJgeLcDhMe3s7ZrN5SJX+FT1Hgr3oUjQapba2lgMHDuByuWhvbycUCtHa2moMsYfDYbZv386aNWt45513jqhfUFVVxfTp0znvvPOYOXMmTU1NZGdn09bWht/vx2q1pgT7jo4OIwHucPn5+UYgdzgcxvK8ZCNGjJBfeKJLgUCAnTt3orUmGAxSU1NDZWVlfzfrhCSmxxoaGlBK9coUlhj8JNiLLvl8PmNJW/K+8XV1dSxevJj333+fVatWpczbp6enM2nSJC6++GImTJhAXl4emZmZ2O12YySgo6PDqGTY2dKlpqYmI+D7/X4sFguFhYVGCdpEe0pKSnC5XLS1taGUMgqJyLCs6EpiE6WEaDTa5/tlhEIhGhsbCYVCuFyuE16nHwgEjIqJie2UJdiLw0mwF13yer3G1yaTie3bt/PjH/+YLVu2pDyuoqKC6dOnM2vWLEaPHo3FYsFqtZKWlobX66WtrY22tjby8/NRSuH3+43et9frPaKHnphLPfXUU40Ssp1lvVut1pOuJS+GH7vdTnl5OXv37jWWVvb1DWJTU5OxAuXQoUOMGTPmhKafEv83Ekmxdru9R9sphgYJ9qJLiWx5p9Np1Ir/xz/+gc1m47zzzuO8885j/vz5uN1usrOzjSH6hoYGQqEQbrcbn89n9KLq6+vJzc3l0KFDtLe3Y7fbjZ3DsrOzjcdZrVaKiopkWZvoFYkaDYlRoP7Izk+MksE/y0WfCLvdzimnnEJLSwsWi8VYfipEMvlNKrqU6GlYrVba29vJz8/nz3/+M0VFRUa99USvPBwOY7fb8fl85Obm0traSkNDAzk5OQSDQVpbW41lT5WVlcbQe01NDeFw2OjZK6U45ZRTJNCLXtefa+1dLhd1dXVEo1GysrJOav8CWRYqjkV+m4ouJTb68Hq9hEIh8vLyOOOMM1J2+wqFQlitVlpbW42a762traSlpZGdnU1mZiZut5tQKITJZDKCeGKO0uFw0NzcTHt7u7G2Xn5xiaEuUdI4EomQlpYmN7eiV8mnS3QpIyODwsJCduzYARy5rSfE5tzz8vJS9oHPzMwkGo1it9txu90opY7ai5JeiRiuZImo6CsS7MUxOZ1O3G73ERvLJKuvrzc2ekkE9sQe7ZIZL4QQ/avPiygrpUYrpTYl/WlVSt2hlLpfKbU/6fiCpOfcp5TarpT6XCk1t6/bPNwlkpkSWfSdFR/RWhtz7uXl5YwYMcJI6hNCCNG/+rxnr7X+HDgLQCllBvYDrwA3Ar/SWj+Y/Hil1DjgauA0YASwXClVpbWWeqh9yOVyceqpp1JbW4vdbje2e02WfFMghBBi4OjvYfyLgB1a6z1dBIgvAS9orQPALqXUdmAS8H4ftVHEuVwuHA4HwWAQl8tFQ0ODsWWtx+MhPz8fp9PZz60UQhyPSCRCS0sL4XAYp9MpeQRDVH8H+6uB55O+/65S6npgA/DvWusmoBhYl/SY6vixIyilbgZuBigrK+uVBg93JpOJtLQ00tLScLvdxjr85Op4QojBo7Gxkb179wKx/8dVVVVDdvvf4azfNj5WStmAy4G/xA/9BjiF2BD/AeAXiYd28nTdyTG01o9rrSdqrSfm5eX1bIPFEZIDvwR6IQan5D0oAoFASrEfMXT0W7AH5gMfaa0PAmitD2qtI1rrKPB7YkP1EOvJlyY9rwSo6dOWCiHEEJU89eZwOKRXP0T15zD+NSQN4SulirTWB+LffgX4NP71a8BzSqlfEkvQqwTW92VDhRBiqMrKyjJq6zscjn4pHSx6X78Ee6VUBnAxcEvS4QeUUmcRG6LfnTintf5MKfUisAUIA7dJJr4QQvQMpdQJ77gnBo9+CfZaaz+Qc9ixb3Tx+J8AP+ntdgkhhBBDUX/O2QshhBCiD0iwF0IIIYY4CfZCCCHEECfBXgghhBjiJNgLIYQQQ5wEeyGEEGKIk2AvhBBCDHES7IUQQoghToK9EEIIMcRJsBdCCCGGOAn2QgghxBDXn7veiUEiGAzS0dGB2WwmLS0Ns9nc300SQghxHCTYi6OKRCI0NDRQU1NDJBLbaNDhcFBcXIzD4ejn1gkhhOguGcYXR9XQ0MC+ffuoqanhN7/5DS+++CI+n49t27bh9/v7u3lCCCG6SXr24gihUIhgMMjf/vY3XnjhBVatWkUkEqGgoIArrrgCi8VCS0sLGRkZ/d1UIYQQ3SDBXhg6OjrYtWsXTz31FC+99BLbt28HwGw2c/HFF3PllVca8/Wtra0UFRX1Z3OFEEJ0kwT7Qaijo4O2tjYAMjIySE9P7/LxwWAQn89HOBzGarUSDAYJh8PYbDacTidpaWls3LiRn//857z22mu0t7cDkJeXxxVXXMGXv/xl8vLyUl7TYpGPjhBCDBb98htbKbUb8AIRIKy1nqiUygb+DFQAu4F/0Vo3xR9/H3BT/PHf01ov64dmDwhNTU3s2rULrTUASinKy8vJyck54rHRaJS2tjZ27dpFKBQiLy+P/fv3E41GAcjMzOSFF17gmWeeYf369cbzzj//fK688kouvvhigsFgp+3weDw9/8MJIYToFf3ZPZupta5P+v5e4G2t9f8ope6Nf3+PUmoccDVwGjACWK6UqtJaR/q+yf0rMcyeCPQAWmv27NlDenq6MYceDodpamrC6/Xi9/sJhUI4nU6am5uJRqPU1dXx5ptv8vTTT9PQ0ADEsuznz5/PwoULqaqqorW1FZfLRVNTk5GJn5CVlYXL5eq7H1wIIcRJGUhjsV8CZsS/fgp4F7gnfvwFrXUA2KWU2g5MAt7vhzb2q7a2tpRAn6C1pq2tjYyMDLTW1NbWcvDgQbKysggEAkBs3n3NmjX89a9/ZeXKlUYAHzVqFP/yL//C17/+dQKBAFlZWcZ71NfXk5WVhVKKQCCA2WzGYrGQl5eH1Wrtux9cCCHESemvYK+BN5VSGvid1vpxoEBrfQBAa31AKZUff2wxsC7pudXxY0dQSt0M3AxQVlbWW20fEDIyMkhLS0NrjVIKkym2itLn83Hw4EHjcT6fj9dff51XXnmFHTt2ALH59rlz53L99ddTVVVFfn4+kUgEu91Oc3MzWVlZxvObmpoAsFqtRKNR3G43mZmZffiTCiGEOFn9Feynaq1r4gH9LaXUP7p4rOrk2JHdWyB+0/A4wMSJEzt9zGCWmZmJUgqHw0EkEqG5uRm3243dbicYDKbMx2/bto3XX3+dl19+2Ui4y8/P56qrruLaa6/F7XYTiURIS0ujsbERk8lETk4OLS0tNDc3k5eXR319vdHLD4VCZGRkUFRUhFKd/ZMIIYQYqPol2Guta+J/1ymlXiE2LH9QKVUU79UXAXXxh1cDpUlPLwFq+rTBA0RaWhplZWV4vV58Ph+5ublEIhFaWlowmUwopVi6dClPPPEEf//7343nTZ48ma985Sv8y7/8Cy0tLTidTlpbW8nIyMBsNhOJRIhEIoRCIeCflfM8Hg9KKWP0wOFwkJaW1l8/vhBCiBPU58FeKZUJmLTW3vjXc4AfA68BC4H/if/9avwprwHPKaV+SSxBrxJYf8QLD1GhUMjoSdfV1dHc3AzEsuGbmprIzMyksbGRZ555hr/85S8pCXcLFizgpptuMobp/X4/drsdn89HRkaGMTSf0NLSgsvlorW1lWg0agzhJ+Tn5yOEEGLw6Y+efQHwSjyAWYDntNZvKKU+BF5USt0E7AWuBNBaf6aUehHYAoSB2/orEz8SidDa2orWGqfT2WtJaoFAAJ/PR2trKy0tLWRlZWGz2Thw4AA2mw2r1YrWmo8//pinn36ad9991wjap5xyCt/61reYNm0adrsdp9NpDOP7/X5ycnLw+XzYbDaAlJ56OBzGZDKRlpZGR0eHcTyxvE/m6oUQYnBSnWV3DwUTJ07UGzZs6NHXrKmp4cCBA0Bs+Vl5eTlms5lQKERzczPhcBiHw4HT6ezydYLBIE1NTfh8PtLT08nKyjIK4/j9fg4cOEB7e7sReNPS0mhubqa9vR2v18uKFSt45pln2LlzJxBLuJszZw5XXnklp59+OllZWcYIgFLKGAUAsNvtxmiAyWQiNzcXn8+XUuve4XBgs9nQWpORkYHL5ZLSuEIIMcAppTZqrSd2dm4gLb0b0KLRqDFEDrEs9eLiYsxmMwcPHjQy4E0mE6NHjz5qcAyHw+zbt88Ixs3NzTQ2NlJVVYXNZqOpqQmTyUQgECAnJ4dwOExraytbtmzhueeeY+nSpUavOz8/n+uuu44rrrgCj8eD0+mkvr4ev9+Pw+HA5/OhtUZrjc1mIxgMEggEiEQiZGdnYzabsdvtuFwuamtr8fl8QCyD32q1UlFRIevphRBiCJBg300mkwm3282hQ4eAWO/XbDajtaalpcV4XDQapaOjg0gkgtfrRSmF0+kkEAiglKKjo8MI9PDPsrNtbW3YbDZaW1tTet+rVq3iySefTEm4O/fcc7n66qu54oorjHX0NpuNUChETk4OTU1NuFwuI3g3NzfjcrlwOBwEAgEsFgtOpxO3221MRTgcDtra2ggGg1gsFjIyMoyhfiGEEIObBPvjUFRURHp6urHePBGoXS6X0dtOrHnftm0bWms8Hg+HDh0iFAqRnZ1tvJbJZCI7O9voaQcCAYLBIBkZGezbt49nn32WF198kcbGRiAWjC+//HKuuOIKKioqgNh0QGFhIYFAwFgiZzab8Xg8RKNRnE4njY2NtLS04PV6yc/Pp6ioqNOMerPZLL14IYQYoiTYHwer1XrEhjAAhYWF2O12wuEwGRkZRhIfxIJ6Ykmb1toYYrfb7dTX12OxWHC5XPj9ft555x2eeuopli9fbiTcVVVVcfXVVzN79mzy8/ONuXSbzUZBQcExkwQ9Ho+R0S89dSGEGJ4k2PcAq9VKfn4+Wmu2b9+esiNccgJkYi7ebrejlDL+PPLII/ztb38zEu6sViuXXHIJ3/rWtxg1ahQZGRmEQiFjWD4nJ4fc3NxurQYwmUzY7fYe/omFEEIMJhLsT5DWGp/PRzQaJTMzE4vFYqxlt1qteDwempubjf3fAbxeL7m5uTQ3N7Nt2zYWLVrEq6++akwBFBQU8JWvfIUvf/nL5ObmopQiMzPTqHZXUFCAzWYjPT1dqtgJIYToNgn2J6ixsZHdu3cDsWV4ZWVlmM1msrKyjDn6rKws3G43WVlZeL1eAoEAr7322hFbyk6ZMoWvf/3rnHPOOUeMClit1pS5fiGEEOJ4SbA/QXV1sWq+iU1jEsVvEglwiVryLpeL3bt38+tf/5pnnnnGWO/udDr5+te/zvz58xk9ejR+vx+Px5OyvE8CvRBCiJ4gwf4EJZameb1eY/94rTW5ubnk5+cTjUZZsmQJDz30EG+//baRcFdZWcmVV17J/PnzOffccwkEAjQ2NuL1emltbSU7OxutNSaTifz8fClmI4QQ4qRJsD9BhYWFNDQ0pKyZ9/l8dHR08Pzzz/P73/+ebdu2AbEe+pw5c/ja177GmWeembIlrcvlwmKxEA6HaWlpobGxEbPZbCTmCSGEECdLgv0JSlSeS5TP3bNnDw888ACvvvqqUegmkXB3xRVXUFlZmTJEX1RUZGTJZ2RkMGrUKNra2ohGo6SlpUkGvRBCiB4jwf4kWK1WNm7cyO9+9zs++ugj4/ikSZO47rrrWLBggVFdz+v1UlxcjFLK2KAmmclkOmZNfSGEEOJESLA/AXv37uW3v/0tf/jDH4zyuS6Xi0suuYSvfvWrVFRUYLfbKS4uxm63G3PxOTk5vbZTnhBCCHE0Euy7KRqNsnz5ch599FFef/11I+HuzDPP5LbbbuOqq66ipaWFuro6TCYTZWVlpKenU1paitZa1sULIYToNxLsj8Hr9fKHP/yB3/zmNykJd1dffTXf+c53mDJlihHIHQ4HeXl5mEymlNK0EuiFEEL0Jwn2x3Do0CHuueceQqEQJSUl3Hrrrdx0000UFBQc8djE/vNCCCHEQCLBvguJ3ehuueUWKioqmDt3LuPGjQMgEAhgNptTKt4JIYQQA5Gpr99QKVWqlFqhlNqqlPpMKXV7/Pj9Sqn9SqlN8T8Lkp5zn1Jqu1Lqc6XU3L5sbzAY5IYbbmDGjBlEo1FCoRB79+7l008/Zfv27bS3t/dlc4QQQojj1h/d0jDw71rrj5RSTmCjUuqt+Llfaa0fTH6wUmoccDVwGjACWK6UqtJaR3q7oWazmdLSUnbt2oVSipKSEtrb24318m1tbbS0tJCent7bTRFCCCFOWJ8He631AeBA/GuvUmorUNzFU74EvKC1DgC7lFLbgUnA+73eWCA7O9uoZJeWlkZra2vK+UQlPCGEEGKg6tdIpZSqAM4GPogf+q5S6hOl1J+UUlnxY8XAvqSnVXOUmwOl1M1KqQ1KqQ2J9e89IS0tzUi8czgclJSUYLfbyc3NxePx9Nj7CCGEEL2h34K9UsoBvATcobVuBX4DnAKcRazn/4vEQzt5uu7sNbXWj2utJ2qtJ+bl5fV8o4n15AsKChg3bhzl5eUpS+yEEEKIgahfgr1Sykos0D+rtX4ZQGt9UGsd0VpHgd8TG6qHWE++NOnpJUBNX7a3MzJ8L4QQYrDoj2x8BfwR2Kq1/mXS8aKkh30F+DT+9WvA1Uopu1JqJFAJrO+r9gohhBCDXX9k408FvgFsVkptih/7f8A1SqmziA3R7wZuAdBaf6aUehHYQiyT/7a+yMQXQgghhor+yMZfTefz8Eu6eM5PgJ/0WqOEEEKIIUwmnoUQQoghToK9EEIIMcRJsBdCCCGGOKV1p0vWBz2l1CFgzzEelgvU90FzhhK5ZsdPrtmJket2/OSaHb+hdM3KtdadFpkZssG+O5RSG7TWE/u7HYOJXLPjJ9fsxMh1O35yzY7fcLlmMowvhBBCDHES7IUQQoghbrgH+8f7uwGDkFyz4yfX7MTIdTt+cs2O37C4ZsN6zl4IIYQYDoZ7z14IIYQY8oZtsFdKzVNKfa6U2q6Uure/29OflFK7lVKblVKblFIb4seylVJvKaW2xf/OSnr8ffHr9rlSam7S8Qnx19mulHoovunRkKGU+pNSqk4p9WnSsR67TvHNnv4cP/6BUqqiT3/AXnCUa3a/Ump//PO2SSm1IOmcXDOlSpVSK5RSW5VSnymlbo8fl8/aUXRxzeSzlqC1HnZ/ADOwAxgF2IC/A+P6u139eD12A7mHHXsAuDf+9b3Az+Jfj4tfLzswMn4dzfFz64HJxPY+WArM7++frYev0zTgHODT3rhOwHeA38a/vhr4c3//zL10ze4H7urksXLNYj9HEXBO/Gsn8EX82shn7fivmXzW4n+Ga89+ErBda71Tax0EXgC+1M9tGmi+BDwV//op4MtJx1/QWge01ruA7cAkFdui2KW1fl/H/jcsSnrOkKC1XgU0Hna4J69T8mv9FbhosI+OHOWaHY1cM0BrfUBr/VH8ay+wFShGPmtH1cU1O5phd82Ga7AvBvYlfV9N1x+MoU4DbyqlNiqlbo4fK9BaH4DYfyQgP378aNeuOP714ceHup68TsZztNZhoAXI6bWW96/vKqU+iQ/zJ4aj5ZodJj5UfDbwAfJZ65bDrhnIZw0YvsG+s7ux4bwsYarW+hxgPnCbUmpaF4892rWTa5rqRK7TcLmGvwFOAc4CDgC/iB+Xa5ZEKeUAXgLu0Fq3dvXQTo4Ny+vWyTWTz1rccA321UBp0vclQE0/taXfaa1r4n/XAa8Qm+Y4GB/SIv53XfzhR7t21fGvDz8+1PXkdTKeo5SyAG66PwQ+aGitD2qtI1rrKPB7Yp83kGtmUEpZiQWtZ7XWL8cPy2etC51dM/ms/dNwDfYfApVKqZFKKRuxZIvX+rlN/UIplamUcia+BuYAnxK7HgvjD1sIvBr/+jXg6nhm6kigElgfH1b0KqXOj89jXZ/0nKGsJ69T8mt9DXgnPm84pCQCVtxXiH3eQK4ZAPGf8Y/AVq31L5NOyWftKI52zeSzlqS/MwT76w+wgFjG5g7gB/3dnn68DqOIZaX+HfgscS2IzUW9DWyL/52d9JwfxK/b5yRl3AMTif1n2gE8Qrxo01D5AzxPbCgwROwu/6aevE5AGvAXYslC64FR/f0z99I1exrYDHxC7BdokVyzlGt2AbHh4U+ATfE/C+SzdkLXTD5r8T9SQU8IIYQY4obrML4QQggxbEiwF0IIIYY4CfZCCCHEECfBXgghhBjiJNgLIYQQQ5wEeyEESqkKlbQzXTce/6RS6mu92SYhRM+RYC+EEEIMcRLshRAJFqXUU/FNQ/6qlMpQSv1QKfWhUupTpdTjne3ydbTHKKXeVUr9TCm1Xin1hVLqwvhxs1Lqwfie4Z8opf4tfnyCUmplfEOmZYdVPxNCnAQJ9kKIhNHA41rr8UArsf27H9Fan6u1Ph1IBy7t5HldPcaitZ4E3AH8KH7sZmJ7iJ8df69n43XNHwa+prWeAPwJ+EmP/4RCDFOW/m6AEGLA2Ke1XhP/+hnge8AupdR/ABlANrGSyv932PNmdvGYxCYuG4GK+Nezgd/q2DahaK0blVKnA6cDb8UHBszEyuwKIXqABHshRMLhtbM18BgwUWu9Tyl1P7H64AalVNoxHhOI/x3hn79vVCfvpYDPtNaTT/aHEEIcSYbxhRAJZUqpRLC9Blgd/7o+vk94Z9n3ad14zOHeBL4d3yYUpVQ2sc1I8hLvr5SyKqVOO8GfQwhxGOnZCyEStgILlVK/I7az2m+ALGK7hu0mtjV0Cq11s1Lq9109phN/AKqAT5RSIeD3WutH4kv5HlJKuYn9bvo1sSkBIcRJkl3vhBBCiCFOhvGFEEKIIU6CvRBCCDHESbAXQgghhjgJ9kIIIcQQJ8FeCCGEGOIk2AshhBBDnAR7IYQQYoiTYC+EEEIMcf8/jBUOOmvHuO8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt_df = df_dd.round({\"balance\": -2}).assign(size=1).groupby(\"balance\").agg({\"pv\":\"mean\", \"size\": \"sum\"}).reset_index()\n",
    "\n",
    "plt.figure(figsize=(8,4))\n",
    "sns.scatterplot(data=plt_df, y=\"pv\", x=\"balance\", size=\"size\", color=\"C5\")\n",
    "plt.plot(plt_df.query(\"balance<5000\")[\"balance\"], m.predict(plt_df.query(\"balance<5000\").assign(balance = lambda d: d[\"balance\"] - 5000)), color=\"C0\", lw=2, label=\"regression\")\n",
    "plt.plot(plt_df.query(\"balance>5000\")[\"balance\"], m.predict(plt_df.query(\"balance>5000\").assign(balance = lambda d: d[\"balance\"] - 5000)), color=\"C0\", lw=2)\n",
    "plt.legend(fontsize=14)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ae47834",
   "metadata": {},
   "source": [
    "### The IV Estimate\n",
    "\n",
    "Errata: The book does not center the discontinuity at zero. Therefore, the intercept for the ITTE regression cannot be interpreted as the ITTE. Here, I'm centering the discontinuity at zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "296983f5",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:14.675146Z",
     "start_time": "2024-01-23T12:22:14.637892Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "732.8534752298891"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def rdd_iv(data, y, t, r, cutoff):\n",
    "    \n",
    "    centered_df = data.assign(**{r: data[r]-cutoff})\n",
    "    \n",
    "    compliance = smf.ols(f\"{t}~{r}*I({r}>0)\", centered_df).fit()\n",
    "    itte = smf.ols(f\"{y}~{r}*I({r}>0)\", centered_df).fit()\n",
    "    \n",
    "    param = f\"I({r} > 0)[T.True]\"\n",
    "    return itte.params[param]/compliance.params[param]\n",
    "\n",
    "\n",
    "rdd_iv(df_dd, y=\"pv\", t=\"prime_card\", r=\"balance\", cutoff=5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "63c00714",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:14.692077Z",
     "start_time": "2024-01-23T12:22:14.684009Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "700.0"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(df_dd\n",
    " .round({\"balance\":-2}) # round to nearest hundred\n",
    " .query(\"balance==5000 & categ=='complier'\")[\"tau\"].mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "8ab77e2c",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:17.650364Z",
     "start_time": "2024-01-23T12:22:14.693750Z"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([655.08214249, 807.83207567])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from joblib import Parallel, delayed\n",
    "from toolz import partial\n",
    "\n",
    "def bootstrap(data, est_fn, rounds=200, seed=123, pcts=[2.5, 97.5]):\n",
    "    np.random.seed(seed)\n",
    "    \n",
    "    stats = Parallel(n_jobs=4)(\n",
    "        delayed(est_fn)(data.sample(frac=1, replace=True))\n",
    "        for _ in range(rounds)\n",
    "    )\n",
    "    \n",
    "    return np.percentile(stats, pcts)\n",
    "\n",
    "\n",
    "bootstrap(df_dd,\n",
    "          partial(rdd_iv, y=\"pv\", t=\"prime_card\", r=\"balance\", cutoff=5000))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32937634",
   "metadata": {},
   "source": [
    "### Bunching\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "a5e34d64",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:17.784941Z",
     "start_time": "2024-01-23T12:22:17.652054Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Balance (rounded to nearest 1000)')"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmoAAAEmCAYAAADSoLz1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAptElEQVR4nO3de5hldX3n+/dHUMQgCtIi0iCMQlSMIrYtaqKoUdHEgDleyHECepxhjjHekkkiM5l4TURNjOMYmOAlgkdF4g00KiLexohCc5GrCBGUFoQ2KGJGUdrv+WP9atiUVUV171pda1W9X8+zn733b+/12d9eVb3Xt9Y1VYUkSZKG507LXYAkSZLmZqMmSZI0UDZqkiRJA2WjJkmSNFA2apIkSQNloyZJkjRQ2y93AX3Zbbfdap999lnuMiRJku7Queee+/2qWjN7fMU2avvssw8bNmxY7jIkSZLuUJJvzzXupk9JkqSBslGTJEkaKBs1SZKkgbJRkyRJGigbNUmSpIGyUZMkSRooGzVJkqSBslGTJEkaqBV7wlutPEkW/d6q6rESSZK2jd7WqCW5a5Kzk3w9ySVJXtPGd01yRpIr2v0uE9Mck+TKJJcneerE+COSXNRee1u2ZIktSZI0Un1u+rwFeGJVPQw4EDg0ycHAK4Ezq2o/4Mz2nCQPBo4ADgAOBY5Lsl3LOh44Gtiv3Q7tsW5JkqRB6K1Rq86P29M7t1sBhwEntvETgcPb48OAk6vqlqq6CrgSWJ9kD2Dnqjqruu1ZJ01MI0mStGL1ejBBku2SXADcAJxRVV8Ddq+q6wDa/b3b2/cErpmYfGMb27M9nj0uSZK0ovXaqFXV5qo6EFhLt3bsIQu8fa79zmqB8V8OSI5OsiHJhk2bNm1xvZIkSUOyTU7PUVU/BL5At2/Z9W1zJu3+hva2jcBeE5OtBa5t42vnGJ/rc06oqnVVtW7NmjVL+U9YcZIs+iZJkpZHn0d9rklyz/Z4R+A3gW8ApwFHtbcdBZzaHp8GHJFkhyT70h00cHbbPHpzkoPb0Z5HTkwjSZK0YvV5HrU9gBPbkZt3Ak6pqk8kOQs4JckLge8AzwaoqkuSnAJcCtwKvLiqNresFwHvAXYEPtVukiRJK1pW6olB161bVxs2bFjuMgZrjCePHWPNkiQtRpJzq2rd7HEvISVJkjRQNmqSJEkDZaMmSZI0UDZqkiRJA2WjJkmSNFA2apIkSQNloyZJkjRQNmqSJEkDZaMmSZI0UDZqkiRJA2WjJkmSNFA2apIkSQNloyZJkjRQNmqSJEkDZaMmSZI0UDZqkiRJA2WjJkmSNFA2apIkSQO1/XIXoJUnyaLfW1U9ViJJ0ri5Rk2SJGmgbNQkSZIGykZNkiRpoGzUJEmSBspGTZIkaaBs1CRJkgaqt0YtyV5JPp/ksiSXJHlZG391ku8muaDdnj4xzTFJrkxyeZKnTow/IslF7bW3ZUvO/yBJkjRSfZ5H7Vbgj6vqvCR3B85NckZ77W+r6q8n35zkwcARwAHAfYHPJtm/qjYDxwNHA18FPgkcCnyqx9olSZKWXW9r1Krquqo6rz2+GbgM2HOBSQ4DTq6qW6rqKuBKYH2SPYCdq+qs6s6OehJweF91S5IkDcU22UctyT7Aw4GvtaE/THJhkncn2aWN7QlcMzHZxja2Z3s8e1ySJGlF671RS7IT8GHg5VX1I7rNmPcHDgSuA/5m5q1zTF4LjM/1WUcn2ZBkw6ZNm6YtXZIkaVn12qgluTNdk/a+qvoIQFVdX1Wbq+oXwDuA9e3tG4G9JiZfC1zbxtfOMf5LquqEqlpXVevWrFmztP8YSZKkbazPoz4DvAu4rKreMjG+x8Tbnglc3B6fBhyRZIck+wL7AWdX1XXAzUkObplHAqf2VbckSdJQ9HnU52OB3wcuSnJBG/svwO8lOZBu8+XVwH8CqKpLkpwCXEp3xOiL2xGfAC8C3gPsSHe0p0d8SpKkFS/dgZQrz7p162rDhg3LXcZgbcmp6Lb0d6Sv7D5rliRpOSU5t6rWzR73ygSSJEkDZaMmSZI0UDZqkiRJA2WjJkmSNFA2apIkSQNloyZJkjRQNmqSJEkDZaMmSZI0UDZqkiRJA2WjJkmSNFA2apIkSQNloyZJkjRQNmqSJEkDZaMmSZI0UDZqkiRJA2WjJkmSNFA2apIkSQNloyZJkjRQNmqSJEkDZaMmSZI0UNsvdwHSckuy6PdWVY+VSJJ0e65RkyRJGijXqEk9cm2dJGkarlGTJEkaKBs1SZKkgdqqRi3JJxbxnr2SfD7JZUkuSfKyNr5rkjOSXNHud5mY5pgkVya5PMlTJ8YfkeSi9trbsiXbkyRJkkZqa9eo/cdFvOdW4I+r6kHAwcCLkzwYeCVwZlXtB5zZntNeOwI4ADgUOC7Jdi3reOBoYL92O3Qr65YkSRqNRTVqSXZM8qszz6vqujuapqquq6rz2uObgcuAPYHDgBPb204EDm+PDwNOrqpbquoq4EpgfZI9gJ2r6qzq9rY+aWIaSZKkFesOG7UkzwAuAD7dnh+Y5LQt+ZAk+wAPB74G7D7T6LX7e7e37QlcMzHZxja2Z3s8e1ySJGlFW8watVcD64EfAlTVBcA+i/2AJDsBHwZeXlU/Wuitc4zVAuNzfdbRSTYk2bBp06bFlihJkjRIi2nUbq2qm7YmPMmd6Zq091XVR9rw9W1zJu3+hja+EdhrYvK1wLVtfO0c47+kqk6oqnVVtW7NmjVbU7IkSdJgLKZRuzjJ/w1sl2S/JP8D+ModTdSOzHwXcFlVvWXipdOAo9rjo4BTJ8aPSLJDkn3pDho4u20evTnJwS3zyIlpJEmSVqzFNGovoTsS8xbg/cBNwMsXMd1jgd8HnpjkgnZ7OnAs8OQkVwBPbs+pqkuAU4BL6faHe3FVbW5ZLwLeSXeAwb8An1rUv06SJGnEckeXrUnyG8BXJpomkhw0c0TnUK1bt642bNiw3GUMVp+XNuore2y5fWdLklaOJOdW1brZ44tZo3Y68Lkku0+MvXPJKpMkSdKcFtOoXQ68GfhCkse0Ma8MIEmS1LPtF/GeqqpPJLkc+GCSdzPP6TEkSZK0dBazRi0AVXUF8BvA44CH9lmUJEmSFrFGraoePvH434DnJNm716okSZI0f6OW5E+r6k1J3jbPW17aU02SJEli4TVql7X7c7dFIZIkSbq9eRu1qvp4uz9xZizJLsAPyxM+SZIk9W7egwmS/EWSB7bHOyT5HN1VAa5P8pvbqkBJkqTVaqGjPp9Ldw416K7JGWAN8Hjgr3quS5IkadVbqFH72cQmzqcCJ1fV5qq6jMWdf02SJElTWKhRuyXJQ5KsAZ4AfGbitbv1W5YkSZIWWjP2MuBDdJs7/7aqrgJI8nTg/G1QmyRJ0qq20FGfXwMeOMf4J4FP9lmUbpMs/rKqHowrSdLKsphLSEmSJGkZ2KhJkiQNlI2aJEnSQN1ho5bkbkn+W5J3tOf7Jfnt/kuTJEla3RazRu0fgFuAR7fnG4HX91aRJEmSgMU1avevqjcBPweoqp/QXaVAkiRJPVpMo/azJDsCBZDk/nRr2CRJktSjxVwK6lXAp4G9krwPeCzw/D6LkiRJ0iIatao6I8l5wMF0mzxfVlXf770ySZKkVW7eRi3JQbOGrmv3eyfZu6rO668sSZIkLbRG7W8WeK2AJy5xLZIkSZqw0LU+nzBNcJJ3A78N3FBVD2ljrwb+I7Cpve2/tGuHkuQY4IXAZuClVXV6G38E8B5gR7prjL6svKilJElaBRZzwtu7JvmjJB9J8uEkL09y10Vkvwc4dI7xv62qA9ttpkl7MHAEcECb5rgk27X3Hw8cDezXbnNlSpIkrTiLOT3HSXQN1P8A3g48GHjvHU1UVV8CblxkHYcBJ1fVLVV1FXAlsD7JHsDOVXVWW4t2EnD4IjMlSZJGbTGn5/jVqnrYxPPPJ/n6FJ/5h0mOBDYAf1xVPwD2BL468Z6Nbezn7fHs8TklOZpu7Rt77733FCVKkiQtv8WsUTs/ycEzT5I8Cvjnrfy844H7AwfSHUU6c8DCXFc6qAXG51RVJ1TVuqpat2bNmq0sUZIkaRgWs0btUcCRSb7Tnu8NXJbkIqCq6qGL/bCqun7mcbvI+yfa043AXhNvXQtc28bXzjEuSZK04i2mUVuynfeT7FFVM+djeyZwcXt8GvD+JG8B7kt30MDZVbU5yc1tjd7XgCPp9pWTJEla8RZzZYJvJ9mFbo3X9hPjC57wNskHgEOA3ZJspLsU1SFJDqTbfHk18J9a1iVJTgEuBW4FXlxVm1vUi7jt9ByfajdJkqQVL3d0SrIkr6O7tue/cNv+YVVVgz7h7bp162rDhg3LXcbUkrl205vblpxerq/cPrPHltt3tiRp5UhyblWtmz2+mE2fzwHuX1U/W/qyJEmSNJ/FHPV5MXDPnuuQJEnSLItZo/YGulN0XAzcMjNYVb/TW1WSJElaVKN2IvBG4CLgF/2WI0mSpBmLadS+X1Vv670SSZIk3c5iGrVzk7yB7lxnk5s+Fzw9hyRJkqazmEbt4e3+4ImxAgZ9eg5JkqSxW8wJb5+wLQqRtHien02SVofFrFEjyW8BBwB3nRmrqtf2VZQkSZIWcR61JP8TeC7wEiDAs4H79VyXJEnSqreYE94+pqqOBH5QVa8BHk133U9JkiT1aDGN2k/a/f9Ocl/g58C+/ZUkSZIkWNw+ap9Ick/gzcB5dEd8vqPPoiRJkrS4oz5f1x5+OMkngLtW1U39liVJkqR5N30meWSS+0w8PxI4BXhdkl23RXGSJEmr2UL7qP098DOAJI8DjgVOAm4CTui/NEmSpNVtoU2f21XVje3xc4ETqurDdJtAL+i9MkmSpFVuoTVq2yWZaeSeBHxu4rVFnShXkiRJW2+hhusDwBeTfJ/uFB3/CyDJA+g2f0qSJKlH8zZqVfWXSc4E9gA+U7ddMPBOdFcpkCRJUo8W3IRZVV+dY+yb/ZUjSZKkGYu5MoEkSZKWgY2aJEnSQNmoSZIkDVRvjVqSdye5IcnFE2O7JjkjyRXtfpeJ145JcmWSy5M8dWL8EUkuaq+9LUn6qlmSJGlI+lyj9h7g0FljrwTOrKr9gDPbc5I8GDgCOKBNc1yS7do0xwNHA/u12+xMSUskyaJvkqT+9daoVdWXgBtnDR8GnNgenwgcPjF+clXdUlVXAVcC65PsAexcVWe104OcNDGNJEnSirat91HbvaquA2j3927jewLXTLxvYxvbsz2ePS5JkrTiDeVggrm2o9QC43OHJEcn2ZBkw6ZNm5asOEmSpOWwrRu169vmTNr9DW18I7DXxPvWAte28bVzjM+pqk6oqnVVtW7NmjVLWrgkSdK2tq0btdOAo9rjo4BTJ8aPSLJDkn3pDho4u20evTnJwe1ozyMnppEkSVrRFryE1DSSfAA4BNgtyUbgVcCxwClJXgh8B3g2QFVdkuQU4FLgVuDFVbW5Rb2I7gjSHYFPtZskSdKKl9uutb6yrFu3rjZs2LDNPm9LTlewJfN8bLl9Zo8tt8/sseVKkhaW5NyqWjd7fCgHE0iSJGkWGzVJkqSBslGTJEkaKBs1SZKkgbJRkyRJGigbNUmSpIGyUZMkSRooGzVJkqSBslGTJEkaKBs1SZKkgbJRkyRJGigbNUmSpIGyUZMkSRooGzVJkqSBslGTJEkaKBs1SZKkgbJRkyRJGigbNUmSpIGyUZMkSRooGzVJkqSBslGTJEkaKBs1SZKkgbJRkyRJGqjtl7sASStfki16f1X1VIkkjYuNmqRR25Im0AZQ0tgsy6bPJFcnuSjJBUk2tLFdk5yR5Ip2v8vE+49JcmWSy5M8dTlqliRJ2taWcx+1J1TVgVW1rj1/JXBmVe0HnNmek+TBwBHAAcChwHFJtluOgiVJkralIR1McBhwYnt8InD4xPjJVXVLVV0FXAms3/blSZIkbVvL1agV8Jkk5yY5uo3tXlXXAbT7e7fxPYFrJqbd2MYkSZJWtOU6mOCxVXVtknsDZyT5xgLvnWtP4Tn3CG5N39EAe++99/RVSpIkLaNlWaNWVde2+xuAj9Jtyrw+yR4A7f6G9vaNwF4Tk68Frp0n94SqWldV69asWdNX+ZIkSdvENm/UkvxKkrvPPAaeAlwMnAYc1d52FHBqe3wacESSHZLsC+wHnL1tq5YkSdr2lmPT5+7AR9u5j7YH3l9Vn05yDnBKkhcC3wGeDVBVlyQ5BbgUuBV4cVVtXoa6JUmStqlt3qhV1beAh80x/q/Ak+aZ5i+Bv+y5NEmSpEEZ0uk5JEmSNMFGTZIkaaBs1CRJkgbKRk2SJGmgbNQkSZIGykZNkiRpoGzUJEmSBspGTZIkaaBs1CRJkgbKRk2SJGmgluNan5I0eO16xItSVT1WImk1c42aJEnSQNmoSZIkDZSbPiVpG9qSTargZlVptXONmiRJ0kDZqEmSJA2UjZokSdJA2ahJkiQNlI2aJEnSQNmoSZIkDZSNmiRJ0kB5HjVJWiG87JW08tioSZIWZAMoLR83fUqSJA3Uqlqj5l+FkiRpTEazRi3JoUkuT3Jlklcudz2SJEl9G8UatSTbAX8HPBnYCJyT5LSqunR5K5Mkba0+L1DvFhStFGNZo7YeuLKqvlVVPwNOBg5b5pokSatMkkXfpKUwijVqwJ7ANRPPNwKPmv2mJEcDR7enP05y+SLzdwO+PytrK8q849w+s8eWu0TZK2IeL1G286L/3DmznRf95y5R9oqZFwPO7TN7bLlbmn2/uQbH0qjN9dv+S+uqq+oE4IQtDk82VNW6rSlsOXL7zDa3/+yx5faZPbbcPrPHlttn9thy+8weW26f2WPLXarssWz63AjsNfF8LXDtMtUiSZK0TYylUTsH2C/JvknuAhwBnLbMNUmSJPVqFJs+q+rWJH8InA5sB7y7qi5Zwo/Y4s2ly5zbZ7a5/WePLbfP7LHl9pk9ttw+s8eW22f22HL7zB5b7pJkx8OSJUmShmksmz4lSZJWHRs1SZKkgbJRkyRJGigbNUmSpIEaxVGfY5HuNNTr6a6kUHTneju7pjxiY2y5fWZbc/+5kpZPkqcCh3P7/9enVtWnV1Nun9mjy12N3+l9LOCSPAU4DrgC+G4bXgs8APiDqvrMasi15nHnTuSP6ousz2xz+88eW25f2UneCuwPnER3onfo/l8fCVxRVS9bDbljrLnXebHaGrUeF5yXAU+rqqtnje8LfLKqHrQacvvMtub+c1vGWxnbF9nIah5b7hhrHum8+GZV7T/HeIBvVtV+qyG3z+yx5QJQVavqBlwG7DPH+L7AZVPkXgFsP8f4XYArV0uuNY87t2V8c57x0C2EBpU7xprHljvGmkc6Ly4E1s8xvh64aLXkjrHmPufFatxHbXtu+wto0neBO0+R+27gnCQnA9e0sb3oLnf1rlWU22e2NfefC/DTJOur6uxZ448EfjrA3D6zze0/e2y5fWY/Hzg+yd25bTm1F/Cj9tpqye0ze2y5q3LT5zHAc4C5FnCnVNUbpsh+EHAY3T4LofthnVZVl05Z86hy+8y25m2SexBwPDDXF84fVNW5Q8odY81jyx1jzWOcFxP592Hi/3VVfW+avLHm9pk9ptxV16hBvwt7aaUY0xdZ39nm9p89tty+s6UZq7JR29aSvLqqXr3ac/vMtub+cyUtnyTnVdVBqz23z+yh5nrC2wlJXt1T9FSrwVdQbp/Z1tx/LknOG1Nun9nm9p89ttw+s/tqesaW22f2UHNdozYhyTOq6uPLXYckaXVLsitQVfWD1ZzbZ/ZYcm3UllBGdtLGvnL7zLbm/nNnfcYovsi2Rba5/WePLXeps5PsDbwJeBLwQ7p933YGPge8smadN3Gl5o6x5j7nxVaf12PMN+CpdEfsnAac2h4fOmXmW4FP0h09+uvtdkQb+++rJdeax53bsvemOyp6E9352q4Ebmhj+wwtd4w1jy13jDWPdF6cBTwX2G5ibLv2f/urqyV3jDX3Oi+mmXiMN/pbcI7qpI195VrzuHNbxvi+yEZW89hyx1jzSOfFvP93p/y+GFXuGGvuc16sxoMJnl5VT6+qk6vqy+12MvBbwNOnyP1pkvVzjC/JSRtHlNtntjX3nwuwW1V9sKo2zwxU1eb2/+ReA8ztM9vc/rPHlttn9rlJjkvyqCT3bbdHJTkOOH8V5Y6x5t7mxarbRy3JhcB/qFlnlG4LvXdV1a9tZe6oTto4xpNBWnP/uS37ZOBG4ERuf1Loo+gWUM8ZUu4Yax5b7hhrHum8uAvwQuY4zyfd8umW1ZA7xpp7nRersFHzjNLbILfPbGvuN3eUX2Qjq3lsuWOseYzzQprLqmvUZvS0gAvdBVgnj8I7u6acyWPL7TPbmvvPlbR8xnaUuEfib4Pc1fid3scCLslTgOPojgD6bhteCzyAbk3dZ1ZDrjWPO3cif1RfZH1mm9t/9thy+8pO8lZgf+AkbtvisxY4km6H9Jethtwx1tzrvFhtjVqPC87LgKfVrHOlJNkX+GRVPWg15PaZbc3957aMtzK2L7KR1Ty23DHWPNJ58c2q2n+O8dAd6b3fasjtM3tsucCqPD3HZcxxnhtgX+CyKXKvALafY/wuwJWrJdeax53bMsZ4SpFR1Ty23DHWPNJ5cSGwfo7x9cBFqyV3jDX3OS+2Z/XZntv+Apr0XeDOU+S+GzinHQ00eRTQEcC7VlFun9nW3H8utFN/1Kwjo1miU4r0kNtntrn9Z48tt8/s5wPHJ5nrYLfnr6LcPrPHlrsqN30eAzyH7gzSsxdwp1TVG6bIfjDwO8w6CqiqLp2y5lHl9pltzdskd4ynFBlVzWPLHWPNY5wXE/mjOUq8z9w+s8eUu+oaNeh3YS+tFGP6Ius729z+s8eW21d226dpNEeJ95XbZ/bYclfjpk9aQ7akTVmSewDH0B0FtKYN30B3LdFjq+qHqyHXmsedO5Ef4H7c9oWzXZLrl+iLbMlzx1jz2HLHWPPY5kUWONgtSS9HiQ8xd4w19zovVtsatR4XnKcDnwNOnPmrqv219XzgSVX15NWQa83jzm05ozulyNhqHlvuGGse6bwY1VHifeX2mT22XGBVHvV5OvBnwH0mxu4DvBI4Y4rcy7fmtZWWa83jzm3T93VkdC+5Y6x5bLljrHmk82JUR4n3lTvGmvucF6tx0+c+VfXGyYHq1kgcm+QFU+R+O8mf0q3huB4gye50aziuWWjCFZbbZ7Y1958L/R0Z3Vdun9nm9p89ttw+s8d2lPhKORJ/b+C5PeQuybxYjY1aXwu459Ktlftiyyvgerprvz1nBLkA3wM+PmXuXNlDnxd9Zved+4WJn99SzQu/1Jcnd9ALi3myh1zz6H7fquoNSU6lO9jt0dx2sNvzaoqD3Vrux+iuTbrUuUteb5/ZPc+LJc+F1bmP2i50C7jDgHu34ZkF3LFV9YMpsh9It5/CV6vqxxPjh9Z0lxVZD1RVnZPkAOBQutXrn9zazHk+571V9ftLmdlyf4PbTvo3zX4hjwK+UVU3Jbkb3c/xIOAS4K+q6qYpsl8KfLSqpl0bNTv3LsDvAd+tqs8meR7wGLqDWU6oqp9Pkf0A4Jl0C4hbgW8CH5hmPkxk93Xqjwcxx4Wsp81t2WM7DUov86LneTyqmvv62fWdrW0ryb2r6oblrmM+q65RW0iSF1TVP2zltC8FXky378KBwMuq6tT22nlVddBW5r4KeBrd2s8z6BqeLwK/CZxeVX+5lbmnzTH8RLod1Kmq39ma3JZ9dlWtb4//A918+RjwFODjVXXsVuZeAjysqm5NcgLwb8CHgSe18d+douabWt6/AO8H/rGqvr+1eRO576P72e0I3AT8CvDRVnOq6qitzH0p8NvAl4CnAxcAP6Br3P6gqr4wbe3atoa+sJhLkntV1b8udx0rRXo+mnuez/xUVT1tK6fdma7etXQ7zH9g4rXjquoPpqjrPsCrgF8AfwG8BPhd4Bt0y9frtjJ31zmGzwMeTvedfONW5p4HfAR4f1V9a2sy5jXNDm4r7QZ8Z4ppLwJ2ao/3ATbQ/TIBnD9l7nbA3ehOprhzG98RuHCK3POA/w84BHh8u7+uPX78lPPx/InH5wBr2uNfYbpLdFw28fi8Wa9dMG3NwJ3omsl3AZuATwNHAXefIvfCdr893Zrb7drzTPnzu2gi627AF9rjvaf5fWsZ9wCOpftC/Nd2u6yN3XOa7AU+81NTTr8z8AbgvcDvzXrtuCly70N3YtO/A+4FvJruUjGnAHtMkbvrHLergV2AXafIPXTWz/Gdrd73A7tPOY+PBXZrjx8BfItuB+pvT/Od0b6L/hz4d0v8O7UO+Hz7ntuL7g/dH7bvpIdPmb0T8Fq6tfk3te+LrwLPnzK3r4PdDprn9gjguilyP9x+Lw6n2yr1YWCHmZ/rlPPi03TN2Svb7/Cfte+3lwCnTpH7C+CqWbeft/tvTZF7FfDXwHeAs4FXAPddkt/lpQgZ0639wOe6XQTcMkXupbOe79R+0d7CFE0Et296zp/12jS5d2q/SGcAB7axrf4lnZX99bbAuRewYb5/z1bk/iPwgvb4H4B17fH+wDlT1jy78bsz3WaNDwCbpsi9mO6on12Am2kLYeCuTHd02EUTX4i7AOdOfuaU82JUC4uW3csCY4QLi/MmHr8TeD3dub5eAXxsynl80cTjzwOPbI/3n/3/fAtze1nAtayn0e16cA3wrDb+JOCsKbNPpduveS3wR8B/A/YDTqTbDWNrc/s6Snwz3daSz89x+8kUuRfMev5fgX+m++6ftlE7f+Lxdxb63C3M/c/t//WvTf4OLsHv2+T/vd+gO33L99o8Pnqq7GmLG9uNbq3Gge3La/K2D3DtFLmfozU8E2PbAycBm6fI/Rpwt/b4ThPj95j2P0LLWUvXAL199n+GKTKvpvtr+6p2f582vtOU/8HuAbyHbvPk1+gWbN+i2xT8sClrPn+B13acIvcVrcZvAy8FzgTeQddovWqK3JfRNQ4n0K35mmlg1wBfmnJejGph0bIvmPV8SRYYI19YzJ4nW11vm/4btNMP0O2HO/naNGvKe1nA3cHP7vytzW3Tf33W83Pa/Z3o9qPd2tzPAH/KxNpPYHe6PxA+O0XuxcB+87x2zRS5lzGxXGpjR9Gtafz2Us1j4PVL9fvWpp9Z7r2F7jJgU6+kmOt7hm5r2KHAP0yVPW1xY7vRbdb69Xlee/+UP/j7zPPaY6fI3WGe8d0mv+SXYL78FlP8JbjIz7gbsO8S5NwdeBjdmpipNudMZO7f47/7vrQ1BMA9gWcB65cg94CW9cAlrndUC4s2fS8LjBEuLDbSreH5Y7o/EDLx2lZvam/Tv6T9bjyRbhPwW4HHAa8B3jtFbi8LOOAsul0Znk33h9LhbfzxTLEGsGV8ZWY5AjyDbn/hmdem+WNmF+CNdE3xD4Ab2+/2G5luk/izgF+d57XDp8h9E/Cbc4wfClwx5Tx+LW13olnjDwA+NE32RNYz6DZZf28Jsk5eiprmzO4r2Js3b+O8zVpY3DhrYbHLFLm9LCza9L0sMEa4sHjVrNvM/qH3AU5agvxDgA/S7dN5EfBJ4GjmONHnFmT2soCj+2PudOBTwAOB/063j9olwGOmzH4o3abVHwJfpv2hR7dG+6VTZj+Q7mCxnWaNH7oEuU/ahrlPW4KfYe810+3v/ZBBz+NpZ6Q3b95Wz422iXUsuWOpedbCYvD1jnEej2Fe0O0ecTndUfJXA4dNvDbNJvy+cl/SR26f2aOcF9NM7M2bt9V1Y4n2Y9xWuWOseWy5Y6x5qPOCfs8eMJrcMdbc57xYjVcmkLSAJBfO9xLdvmqDyu0z29z+s8eW23P2dtVOll5VVyc5BPhQkvu17NWS22f22HJt1CT9kt2Bp9LtzDwpdDtRDy23z2xz+88eW26f2d9LcmBVXQBQVT9O8tt0l6z6tVWUO8aae5sXNmqSZvsE3Sr8C2a/kOQLA8ztM9vc/rPHlttn9pF0l4P7P6rqVuDIJH+/inL7zB5brpeQkiRJGqo7LXcBkiRJmpuNmiRJ0kDZqEmSJA2UjZq0AiXZnOSCJF9Pcl6Sxyximh9vi9rm+Nwk+VySnbfx516dZLcteP/zk7x9jvFDFjN/hyjJ4UkePM9rj2u/O7cmedas145KckW7HTUxvm+Sr7XxDya5SxtPkrcluTLJhUkOauN3SfKlJB7YJs3DRk1amX5SVQdW1cOAY4A3LHdBC3g63TU1fzQ52BbuY/iOOgQYRKOWZLstnORwYM5GDfgO8Hzg/bM+Y1e6y1M9ClgPvCrJLu3lNwJ/W1X70Z264oVt/GnAfu12NHA8QFX9DDgTeO4W1i2tGmP4EpQ0nZ1p53tKslOSM9uakouSHDb7zfO9J8k+SS5L8o4klyT5TJId22sPSPLZiTV492/jf5LknLYW5TXz1Pc84NRZn3EccB6wV5I3J7m41fLc9r5Dknxioua3J3l+e3x1ktdM1P/ANn6vVvP57XD5TEz/75Oc3dZC/v1Mw5PkBUm+meSLwGPnmFf7AP8v8Io27W8kuV+bfxe2+73nmO7VSd6d5AtJvpXkpYuo5fgkG9q8f83E+69O8hdJvgw8O8lTkpzV/v3/mGSn9r5jk1za6vrrthbwd4A3t8+6/2SNVXV1VV0I/GJW+U8FzqiqG6vqB8AZwKFJQnfB9g+1951I1wgCHEZ3rdGqqq8C90yyR3vtY+13QNIcbNSklWnHtvD9BvBO4HVt/KfAM6vqIOAJwN+0Beykhd6zH/B3VXUA3QWp/682/r42/jC6tUvXJXlKe/964EDgEUkeN0etjwXOnXj+q3QL9YcD69q0D6O7UPWbJxbwC/l+q/944D+3sVcBX265pwF7AyR5EN0ancdW1YHAZuB57XNe0+p7MnOseaqqq4H/SbcW6cCq+l/A21v9D23z5W3z1PhAuqZnZq3Uneerpb3/v1bVOroLgj8+yUMnsn5aVb8OfBb4c7oL1B9EdxmbP2prwZ4JHNDqen1VfaXNhz9ptf/LHc7Vzp7ANRPPN7axewE/bOeOmhxfaBqAi4FHLvKzpVXH/QKkleknbUFPkkcDJyV5CN1apL9qDdMv6BaWuwPfm5h2vvcAXDVxks9zgX2S3B3Ys6o+ClBVP22f+xTgKcD57f070TVuX5pV665VdfPE82+3tS4Avw58oKo2A9e3NVuPBH7Ewj4yUePvtsePm3lcVf+UZOas8k8CHgGc0/rRHYEb6DbtfaGqNrV/zweB/e/gcwEePfGZ7wXeNM/7/qmqbgFuSXID3TyerxaA5yQ5mu57ew+6xnHmUkYfbPcHt/F/btPfBTiLbn79FHhnkn+iO2Hr1prrcji1wPhC01BVm5P8LMndZ/0eSMJGTVrxquqsdDvNr6HbH2wN8Iiq+nmSq4G7zprkeQu855aJ922mayTmu45dgDdU1R2dlfvWJHeqqplNbP82K2POabj9FoHZ/4aZOjdz+++5uc7wHeDEqjrmdoPJ4fO8f0vNlzF7Xm6/QC370q0ZfGRV/SDJe7j9v3lmnoVus+Tvzf6wJOvpGsEjgD+k20y5NTbS7Zc3Yy3wBeD7dJs0t29r1dYC105Ms9esaa6deL4DXSMpaRY3fUorXNtHazvgX4F7ADe0BuwJwP3mmGQx7/k/2kEAG1tjQ5IdktwNOB34fyb2kdozyb3niLgc+HfzxH8JeG6S7ZKsoVsrdjbwbeDB7bPuQdeA3JEv0TYjJnkaMLMD/JnAs2ZqS7Jrugspfw04pO3bdmfg2fPk3gzcfeL5V+iaIdrnfXkRtc2Yr5ad6Zqxm5LsTrdz/ly+Cjw2yQPa9HdLsn/7Gdyjqj4JvJxuc/JctS/G6cBTkuyS7iCCpwCnV3eZm88DM0eIHkXb95BuE+uR6RwM3FRV17Ua7wVsqqqfb2Ed0qrgGjVpZdoxyQXtcYCj2iam9wEfT7IBuAD4xhzTLuY9s/0+8PdJXgv8HHh2VX2m7XN1VtsM92Pg33PbprwZ/0S3hubKOXI/Srcp8et0a6b+tKq+B5DkFLpNf1dw2+bVhbwG+ECS84Av0h3VSFVdmuTPgc+kO8r058CLq+qrSV5Nt+nwOrqDG+Y6qvLjwIfSHXTxEuClwLuT/AmwCXjBImpjEbWcD1wCfAv453mm35TuoIoPJNmhDf85XUN2apK70v0+vKK9djLwjnYww7Mm91NL8ki6+b8L8Iwkr6mqA6rqxiSvA85pb31tVd3YHv8ZcHKS19P9TN7Vxj9Jtzb3SuB/z5onT2ivS5qD1/qUtKzaTvsnVdWTl7sWbXtJPgIcU1WXL3ct0hC56VPSsmqbwN6RbXzCWy2/dCfE/ZhNmjQ/16hJkiQNlGvUJEmSBspGTZIkaaBs1CRJkgbKRk2SJGmgbNQkSZIG6v8HSsp7L2/GF64AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10,4))\n",
    "df_dd.round({\"balance\":-3}).groupby(\"balance\").size().plot.bar()\n",
    "plt.ylabel(\"Sample Size.\")\n",
    "plt.xlabel(\"Balance (rounded to nearest 1000)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38edfa54",
   "metadata": {},
   "source": [
    "## Key Ideas\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "5c420823",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-01-23T12:22:17.958952Z",
     "start_time": "2024-01-23T12:22:17.786257Z"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.50.0 (0)\n",
       " -->\n",
       "<!-- Pages: 1 -->\n",
       "<svg width=\"301pt\" height=\"98pt\"\n",
       " viewBox=\"0.00 0.00 300.98 98.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 94)\">\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-94 296.98,-94 296.98,4 -4,4\"/>\n",
       "<!-- U -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>U</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"37.7\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"37.7\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">U</text>\n",
       "</g>\n",
       "<!-- School -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>School</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"146.49\" cy=\"-64\" rx=\"35.19\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"146.49\" y=\"-60.3\" font-family=\"Times,serif\" font-size=\"14.00\">School</text>\n",
       "</g>\n",
       "<!-- U&#45;&gt;School -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>U&#45;&gt;School</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M60.82,-27.53C74.92,-33.6 93.5,-41.61 109.67,-48.57\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"108.29,-51.79 118.86,-52.53 111.06,-45.36 108.29,-51.79\"/>\n",
       "</g>\n",
       "<!-- Income -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>Income</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"255.28\" cy=\"-41\" rx=\"37.89\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"255.28\" y=\"-37.3\" font-family=\"Times,serif\" font-size=\"14.00\">Income</text>\n",
       "</g>\n",
       "<!-- U&#45;&gt;Income -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>U&#45;&gt;Income</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M64.7,-20.77C100.16,-24.56 164.26,-31.4 208.02,-36.06\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"207.81,-39.56 218.13,-37.14 208.56,-32.6 207.81,-39.56\"/>\n",
       "</g>\n",
       "<!-- School&#45;&gt;Income -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>School&#45;&gt;Income</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M179.24,-57.16C189.12,-55.03 200.19,-52.65 210.73,-50.38\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"211.57,-53.78 220.61,-48.25 210.1,-46.94 211.57,-53.78\"/>\n",
       "</g>\n",
       "<!-- Quarter -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>Quarter</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"37.7\" cy=\"-72\" rx=\"37.89\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"37.7\" y=\"-68.3\" font-family=\"Times,serif\" font-size=\"14.00\">Quarter</text>\n",
       "</g>\n",
       "<!-- Quarter&#45;&gt;School -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>Quarter&#45;&gt;School</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M75.15,-69.27C83.66,-68.63 92.82,-67.95 101.61,-67.29\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"102.01,-70.77 111.72,-66.53 101.48,-63.79 102.01,-70.77\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.graphs.Digraph at 0x7f84bb6f7490>"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gr = Digraph(format=\"png\", graph_attr={\"rankdir\":\"LR\"})\n",
    "\n",
    "gr.edge(\"U\", \"School\")\n",
    "gr.edge(\"U\", \"Income\")\n",
    "gr.edge(\"Quarter\", \"School\")\n",
    "gr.edge(\"School\", \"Income\")\n",
    "gr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3636bfae",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python [conda env:root] *",
   "language": "python",
   "name": "conda-root-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.12"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
